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Preface 


This  thesis  employs  the  application  of  the  methodology  developed  by  Johnson, 
Bauer,  Moore,  and  Grant  to  a  large  scale  linear  programming  (LP)  model,  specifi¬ 
cally  AMC’s  STORM.  The  methodology  of  Johnson,  et  al.  utilizes  response  surface 
methology  and  kriging.  It  basically  accomplishes  an  optimality  analysis  for  LP  mod¬ 
els  by  developing  first  or  second  order  metamodels  which  describe  the  relationships 
between  the  optimal  objective  function  value  and  the  RHS  vectors  of  the  LPs. 

This  research  may  benefit  the  analysts  of  AMC  because  the  developed  meta¬ 
models  can  provide  some  valuable  insights  about  the  system  when  the  resources  of 
the  channel  cargo  system  are  changed. 
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Abstract 

A  methodology  for  optimality  analysis  of  linear  programs  was  developed  by 
Johnson,  Bauer,  Moore,  and  Grant  to  create  metamodels  using  response  surface 
methodology  techniques  such  as  experimental  design  and  least  squares  regression, 
and  a  geostatistical  estimation  technique,  namely  kriging.  Metamodels  have  the  form 
of  a  simple  polynomial,  and  they  predict  the  optimal  objective  function  value  of  an 
LP  for  various  levels  of  the  constraints.  They  eliminate  the  necessity  of  determining 
which  critical  region  contains  the  right-hand-side  (RHS)  vector  of  interest  since  they 
are  valid  over  multiple  critical  regions. 

The  methodology  of  Johnson,  et  al.  can  be  applied  to  large  scale  linear  pro¬ 
gramming  models.  The  developed  metamodels  of  the  large  scale  LP  can  provide  some 
useful  information  about  the  relationships  between  the  objective  function  value  and 
the  RHS  vector  of  interest. 
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RIGHT-HAND-SIDE  MULTIDIMENSIONAL  OPTIMALITY 
ANALYSIS  OF  A  LARGE  SCALE  LINEAR  PROGRAM 
USING  METAMODELLING  TECHNIQUES 


/,  Introduction 

This  research  supports  the  efforts  of  Johnson,  Bauer,  Moore,  and  Grant  (21), 
and  investigates  the  application  of  the  methodology  they  have  developed  to  a  large 
scale  linear  program. 

The  methodology  uses  experimental  design,  linear  regression,  and  kriging  tech¬ 
niques  to  perform  an  optimality  analysis.  This  study  is  sponsored  by  the  United 
States  Air  Force  Air  Mobility  Command  (AMC).  The  first  chapter  provides  the 
background,  the  research  objective,  the  scope  of  the  study,  and  a  summary  of  cur¬ 
rent  knowledge. 

1.1  Background 

Linear  programming  (LP)  is  a  tool  for  solving  optimization  problems.  Since 
the  development  of  the  simplex  algorithm  for  solving  LPs  by  George  Dantzig  in  1947, 
and  later  the  development  of  interior  point  methods,  LP  has  been  used  to  solve  the 
optimization  problems  of  industries  and  government.  In  a  survey  by  Fortune  maga¬ 
zine,  85  percent  of  500  responding  firms  indicated  that  they  used  linear  programming 
for  optimization  purposes  (35:47). 

It  is  difficult  to  give  a  precise  definition  of  a  “large”  LP  model  since  the  defini¬ 
tion  changes  according  to  the  user  and  the  tools  used  for  solving  and  analyzing  these 
large  models.  Brook,  Kendrick,  and  Meeraus  (9:166-167)  describe  large  models  as 
follows: 
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What  is  a  large  model?  The  answer  varies.  It  is  one  that  takes  a  lot  of 
time  (or  money)  to  solve.  It  is  one  that  “just  fits”  into  memory  available 
on  your  machine.  It  is  one  that  has  more  than  a  few  hundred  lines  of 
assignment  and  equations  when  written  in  GAMS.  Briefly,  any  model 
that  is  expensive  to  solve  or  difficult  to  manage,  or  whose  details  are  so 
overwhelming  that  it  hard  to  keep  track  of  them  is  large. 

Approaches  to  solving  large  scale  LPs  were  initiated  with  the  publication  of  the 
Dantzig- Wolfe  decomposition  principle  in  1960.  However,  there  was  significant  work 
in  this  area  even  before  that.  Since  then,  the  growth  of  studies  has  been  explosive 
(26:144).  As  the  size  of  LPs  increased,  the  difficulty  of  optimality  analysis,  as  de¬ 
scribed  below,  has  increased.  Although,  today  we  are  able  to  solve  large  LPs  with 
millions  of  variables  and  hundred  thousands  of  constraints  by  using  highly  developed 
computer  packages,  we  still  need  efficient  optimality  analysis  methods. 

The  right-hand-side  (RHS)  vectors  of  an  LP  may  be  changed  for  a  variety  of 
reasons:  we  may  want  to  conduct  a  sensitivity  analysis  to  check  the  preciseness  of  the 
RHS  vectors,  and  whether  or  not  it  matters  if  the  RHS  vectors  are  perturbed.  Also, 
we  may  want  to  update  the  LP  when  additional  (reduced)  resources  become  available 
(unavailable).  In  this  case,  optimality  analysis  comes  into  play.  Optimality  analysis 
is  performed  to  determine  the  effect  on  the  optimal  solution  when  the  right-hand- 
side  vector  (or  the  objective  function  coefficient)  is  changed.  Optimality  analysis 
generally  involves  multiple  critical  regions  with  different  optimal  bases.  This  differs 
from  “post-optimality  analysis”  and  “sensitivity  analysis”  since  those  analyses  deal 
with  only  one  critical  region  (21).  Optimality  analysis  of  large  scale  LPs  across 
multiple  critical  regions  is  more  difficult  than  situations  dealing  with  only  one  critical 
region.  Thus,  identifying  which  critical  region  contains  a  particular  right-hand-side 
vector  creates  a  burden  for  the  analyst. 

Unlike  other  multiparametric  techniques  developed  by  Gal  and  Nedoma  (14) 
Mulvey,  et  al.  (28)  and  Wagner  (32),  the  methodology  proposed  by  Johnson,  et  al. 
creates  metamodels  by  using  experimental  design,  linear  regression,  and  kriging  to 
perform  optimality  analysis.  This  methodology  eliminates  the  necessity  of  determin- 
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ing  which  critical  region  contains  a  right-hand-side  vector  of  interest.  Furthermore, 
it  insures  efficient  data  requirements,  and  it  allows  us  to  understand  the  key  rela¬ 
tionships  between  the  objective  function  and  the  constraints.  The  method  develops 
metamodels  that  can  provide  valuable  insights  about  the  system.  Sometimes  the 
metamodels  are  able  to  represent  the  whole  system  (21). 

1.2  Research  Objectives 

We  will  employ  a  large  scale  LP  where  the  right-hand-side  (RHS)  vectors  can 
vary  over  a  range  of  values.  A  multidimensional  metamodelling  technique  for  esti¬ 
mating  the  objective  function  over  the  right-hand-side  vector  has  been  developed  by 
Johnson,  et  al.  The  problem  is  to  verify  whether  or  not  this  technique  is  valid  for 
large  scale  linear  programs.  Since  metamodels  are  really  time  and  effort  savers,  the 
analyst  will  be  able  to  observe  the  response  of  the  optimal  objective  function  value 
of  a  particular  LP  very  easily  and  efficiently  when  the  levels  of  the  constraints  of 
this  LP  are  changed. 

1.3  Scope  of  the  Research 

This  research  performs  optimality  analysis  only  on  the  right-hand-side  vec¬ 
tor.  We  intend  to  create  metamodels  with  only  first  order  polynomials  unless  the 
higher  order  polynomials  (such  as  second  order)  are  needed.  We  basically  apply 
2^“*^  fractional  and  2*  full  fractional  designs  to  create  metamodels  by  least  squares 
regression.  Furthermore,  in  this  research  we  apply  only  ordinary  kriging  techniques 
(point  kriging  and  block  kriging)  assuming  that  there  is  no  drift  in  the  data  set. 

For  the  evaluation  of  metamodels,  three  primary  measures  are  used:  mean 
squared  error  (MSE),  percentage  of  predictions  closer  to  the  true  optimal  solution 
(PC),  and  mean  absolute  percentage  error  (MAPE). 
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1.4  Summary  of  the  Current  Knowledge 

In  this  thesis,  STORM  was  chosen  as  a  large  scale  LP  for  the  application  of 
the  methodology  of  Johnson  et  al.  This  section  first  introduces  STORM,  and  then 
gives  a  brief  summary  of  the  approaches  to  sensitivity  analysis. 

1.4.1  Strategic  Transport  Optimal  Routing  Model  (STORM)  .  The  Strate¬ 
gic  Transport  Optimal  Routing  Model  (STORM)  is  based  on  a  model  built  by  Barton 
and  Gumaer  (1967)  for  Lockheed  to  analyze  the  peacetime  employment  of  the  new 
C-5  cargo  plane.  STORM  was  developed  at  the  Air  Mobility  Command  (then  the 
Military  Airlift  Command)  to  assist  in  a  major  study  of  the  entire  scheduled  cargo 
system  that  must  provide  two  main  types  of  service  to  its  overseas  customers.  The 
first  is  to  provide  sufficient  cargo  capacity  for  a  given  period  of  time  (usually  for  one 
month)  to  meet  all  demands  for  cargo  movement  between  the  pairs  of  cities  in  the 
system.  This  cargo  capacity  is  known  as  the  cargo  requirement.  The  second  is  to 
provide  a  minimum  number  of  flights  per  month  between  certain  cities.  This  number 
is  called  the  frequency  requirement.  The  basic  purpose  of  STORM  “is  to  select  the 
mix  of  routes  and  aircraft  that  will  meet  the  monthly  cargo  and  frequency  require¬ 
ments  while  minimizing  the  cost  of  cargo  handling,  military  aircraft  operations,  and 
commercial  aircraft  leasing”  (1:2). 

The  set  of  routes,  maximum  payload,  and  total  flying  hours  for  each  type  of 
plane  are  the  main  resources  available  from  which  STORM  constructs  a  feasible  cargo 
movement  plan.  Routes  are  the  sequences  of  legs  to  be  flown  by  a  single  aircraft. 
Each  aircraft’s  maximum  payload  is  an  average  based  on  the  fuel  loads  and  types 
of  palletized  cargo  that  are  generally  carried  on  the  planned  missions.  The  flying 
hour  limit  for  military  aircraft  is  derived  from  the  Air  Force  flying  program  which  is 
necessary  to  maintain  proficiency  in  worldwide  operations  and  to  train  the  crews. 

Versions  of  STORM  for  UNIX  Workstations  have  been  developed  using  the 
GAMS  modelling  language  (9)  which  makes  data  management  and  programming 
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very  easy.  GAMS  also  allows  the  analyst  to  modify  the  model  quickly  for  specific 
analyses,  investigating  specific  questions,  or  enforcing  operational  considerations  lo¬ 
cally.  For  detailed  information  about  the  modelling  of  STORM  using  GAMS,  refer 
to  Appendix  A. 

In  the  next  sections  we  divide  the  sensitivity  analysis  approaches  into  two  cases 
according  to  the  critical  regions  they  are  dealing  with. 

I.4.S  The  Case  When  the  Candidate  Set  is  In  a  Single  Critical  Region.  If 
we  do  not  want  to  change  the  optimal  basis,  or  we  are  certain  that  potential  changes 
in  a  set  of  RHS  vectors  will  occur  in  a  single  region,  then  we  can  use  one  of  the 
techniques  developed  for  multidimensional  sensitivity  analysis.  Some  examples  of 
these  techniques  are  summarized  below. 

Higher  Dimensional  Cuts.  A  cut  is  a  line  which  divides  the  critical  region  of 
the  candidate  set  into  two  parts.  Making  a  one-dimensional  cut  through  the  critical 
region  of  the  RHS  vector  and  characterizing  the  end  points  of  this  cut  was  proposed 
by  Gass  in  1954.  One-at-a-time  restriction  is  a  major  limitation  of  this  approach 
(33:14).  The  RHS  vectors  often  vary  simultaneously  and  dependently  so  that  one¬ 
dimensional  cuts  are  not  applicable.  Rather  than  making  a  one-dimensional  cut  of 
a  critical  region,  another  approach  to  simultaneous  perturbations  is  to  make  higher 
dimensional  cuts  (see  Gass  and  Saaty  p.26)  (15). 

100  Percent  Rule.  As  proposed  by  Bradley,  Hax  and  Magnanti  (8:7),  one¬ 
dimensional  bounds  can  be  used  indirectly  to  give  an  approximation  to  a  critical 
region  by  using  a  procedure  called  the  100  Percent  Rule.  This  rule  requires  a  spec¬ 
ification  of  the  increase  or  decrease  in  the  direction  for  each  RHS  vector.  This 
guarantees  that  the  same  basis  remains  optimal  as  long  as  the  sum  of  the  fractions 
which  correspond  to  the  percent  of  maximum  change  in  each  direction  is  less  than 
or  equal  to  one. 
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The  Tolerance  approach.  This  approach  was  described  by  Wendell  (34).  In 
contrast  to  “ordinary”  sensitivity  analysis  in  linear  programming,  it  considers  simul¬ 
taneous  and  independent  changes  in  the  objective  function  coefficients  and  in  the 
RHS  vectors.  The  tolerance  approach  yields  a  maximum  tolerance  percentage  such 
that,  as  long  as  selected  coefficients  or  vectors  are  accurate  within  that  percentage 
of  their  estimated  values,  the  same  basis  remains  optimal. 

A  Perturbation  Analysis  Approach  by  Arsham  and  Oblak  (3).  This  technique 
allows  simultaneous  perturbation  of  the  RHS  vector  and/or  the  cost  coefficients 
and/or  the  matrix  coefficients  from  their  nominal  values.  Using  the  perturbed  opti¬ 
mal  solution,  the  current  optimal  solution  is  updated  for  various  points  within  the 
critical  region. 

Convex  Bounds.  The  adaptation  of  the  exploitation  of  convex  bounds  to  es¬ 
timate  the  effects  of  large  perturbations  of  the  RHS  vectors  or  objective  function 
coefficients  in  an  LP  was  derived  from  the  properties  of  convexity  by  Fiacco  (13). 
He  describes  iterative  procedures  for  choosing  additional  points  to  improve  these 
bounds.  This  boundary  approach  helps  to  address  the  computational  difficulties  of 
computing  the  optimal  objective  function  value  for  higher  dimensional  cuts. 

1.4.3  Approaches  Dealing  With  Multiple  Critical  Regions.  If  the  set  of 
RHS  vectors  contains  vectors  from  more  than  one  critical  region,  analysis  of  an 
LP  becomes  more  difficult.  Techniques  developed  for  this  analysis  are  summarized 
below. 

Gal  and  Nedoma’s  Methodology  for  Multiparametric  Programming.  In  their  pa¬ 
per,  Gal  and  Nedoma  (14)  present  an  effective  method  for  finding  all  non-overlapping 
critical  regions  that  cover  a  prespecified  admissable  set.  This  method  uses  an  algo¬ 
rithm  to  determine  which  critical  region  contains  a  particular  RHS  vector.  The 
algorithm  is  developed  so  that  it  uses  an  iterative  approach  and  utilizes  one  of  the 
techniques  described  in  section  1.4.1,  for  example,  the  tolerance  approach  of  Wendell. 
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The  Methodology  of  Mulvey,  Vanderhei,  and  Zenios  (28).  This  methodology 
finds  robust  solutions  to  optimization  problems.  Mulvey,  et  al.  use  various  scenarios 
to  identify  the  uncertainty  in  model  parameters.  They  employ  a  multi-objective  pro¬ 
gramming  framework  to  resolve  conflicting  objectives.  The  solution  to  the  problem 
is  “solution  robust”  if  it  remains  close  to  optimal  across  various  scenarios.  In  con¬ 
trast,  a  solution  is  “model  robust”  if  it  remains  “almost”  feasible  across  the  scenarios 
applied.  Basically  Mulvey,  et  al.  seek  a  simple  robust  solution  in  reaction  to  noisy 
parameters  (21). 

Methodology  of  Wagner  (32).  Wagner  uses  Monte  Carlo  methods  in  his  method¬ 
ology,  and  assumes  knowledge  of  the  joint  probability  distribution  of  a  vector  of 
parameters  of  interest.  He  takes  samples  from  this  distribution,  solves  the  resulting 
optimization  problem,  and  produces  a  data  set  that  serves  as  an  input  to  a  regression 
routine.  Wagner  basically  seeks  a  characterization  of  the  optimal  solution  as  a  func¬ 
tion  of  parameters  which  vary  according  to  a  prespecified  probability  distribution 
(21:15).  The  methodology  proposed  by  Johnson,  Bauer,  Moore,  and  Grant  can  be 
considered  similar  to  Wagner’s  methodology  in  purpose.  However,  the  methodology 
suggested  by  Johnson,  et  al.  uses  response  surface  methods  (specifically  experi¬ 
mental  design  and  least  squares  regression)  and  kriging  instead  of  joint  probability 
distributions  for  Monte  Carlo  methods.  Johnson,  et  al.  develops  metamodels  whose 
form  is  a  polynomial  describing  the  optimal  objective  function  value  as  a  function 
of  the  RHS  vector. 
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11.  Literature  Search  and  Review 


This  chapter  discusses  some  basic  aspects  of  linear  programming  and  the  tech¬ 
niques  used  in  the  methodology  of  Johnson,  et  al.  which  we  will  apply  to  a  large 
scale  system,  specifically  AMC’s  STORM  model.  These  techniques  include  experi¬ 
mental  design,  least  squares  regression,  and  kriging.  Other  approaches  to  optimality 
analysis  were  summarized  in  Section  1.4  of  Chapter  1. 

2.1  Linear  Programming 

The  general  linear  program  is  described  by  Hadley  (17:4)  as  follows: 

Given  a  set  of  m  linear  inequalities  or  equations  in  r  variables,  we  wish 
to  find  non-negative  values  of  these  variables  which  will  satisfy  the  con¬ 
straints  and  maximize  or  minimize  some  linear  function  of  the  variables. 

To  represent  an  optimization  problem  as  a  linear  program,  assumptions  about  pro¬ 
portionality,  additivity,  divisibility  and  deterministic  coefficients  are  made  implicitly 
during  formulation  (4:3-4).  With  the  inclusion  of  slack  variables  and  simple  trans¬ 
formations,  an  LP  can  be  stated  in  matrix  form  as: 

Minimize  cx 

Subject  to  Ax  =  b  (2.1) 

X  >  0 

where  c  is  a  1  x  n  vector  of  objective  function  coefficients,  x  is  an  n  X  1  vector  of 
decision  and  slack  variables,  A  is  an  m  x  n  matrix  of  the  constraint  coefficients,  and 
b  is  an  m  X  1  vector  of  the  right-hand-sides  of  the  constraints. 

2.1.1  Existence  of  Solutions.  For  any  right-hand-side  vector,  b  G  S?*”,  four 
possible  cases  may  arise  for  the  linear  program. 
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1.  Unique  Finite  Optimal  Solution.  When  an  optimal  finite  solution  is  unique, 
it  occurs  at  an  extreme  point.  Figures  2.1a  and  2.1b  show  a  unique  optimal 
solution.  In  Figure  2.1a.,  the  feasible  region  is  bounded.  In  Figure  2.1b,  the 
feasible  region  is  not  bounded. 


Figure  2.1  Unique  finite  optimal  solution:  (a)  Bounded  region,  (b)  Unbounded 
region. 

2.  Alternate  Finite  Optimal  Solutions.  Figure  2.2a  shows  this  case  when  the 
feasible  region  is  bounded.  The  two  extreme  points,  and  are  optimal, 
as  is  any  point  on  the  line  segment  joining  them.  In  Figure  2.2b,  the  feasible 
region  is  unbounded  but  the  value  of  the  objective  function  at  optimality  is 
finite. 


Figure  2.2  Alternate  finite  optima:  (a)  Bounded  region,  (b)  Unbounded  region. 
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3.  Unbounded  Optimal  Solution  Value.  This  case  is  illustrated  in  Figure  ^.5  where 
both  the  feasible  region  and  the  optimal  solution  value  are  unbounded.  Any 
point  on  the  “ray”  with  vertex  x*  is  feasible.  So,  the  value  of  the  solution  is 
unbounded.  No  optimal  solution  exists. 


4.  Empty  Feasible  Region.  In  this  case,  the  system  of  equations  and/or  inequalities 
defining  the  feasible  region  is  inconsistent  (4:18).  In  Figure  2.4,  it  is  clear  that 
there  is  no  point  (a;i,a;2)  satisfying  the  inequalities. 


Figure  2.4  Empty  feasible  region. 
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2.1.2  Interpretation  of  Feasibility.  Consider  the  linear  programming  prob¬ 
lem  given  by  (2.1).  For  the  m  X  n  A  matrix,  the  jth  column  is  denoted  by  aj.  The 

problem  can  be  written  as  follows: 

Minimize  Z)j=i  CjXj 

Subject  to  =  b  (2-2) 

Xj  >  0  ;=l,2,...,n 


Given  the  vectors  ai,  a2, ..,  a,i,  we  are  trying  to  find  nonnegative  scalars  Xi,X2,  ...,Xn 
such  that  CjiCj  is  minimized,  and  SLjXj  =  b.  However,  note  that  the  collec¬ 
tion  of  vectors  of  the  form  X)j_i  SLjXj,  where  Xi,X2, ...,  Xn  >  0,  is  the  cone  generated 
by  ai,a2,...,a„  (4:20).  This  cone  is  also  called  the  positive  cone  of  A,  and  is  de¬ 
noted  by  POS  (A)  (29).  The  problem  has  a  feasible  solution  if  and  only  if  the 
vector  b  G  POS{A)  (see  Figure  2.5a  and  2.5b).  Because  the  vector  b  reflects  some 
requirements  to  be  satisfied,  Bazaraa,  et  al.  call  Figure  2.5  a.  requirement  space 
(4:20). 


Figure  2.5  Illustration  of  the  requirement  space:  (a)  System  is  feasible,  (b)  System 
is  inconsistent. 

The  requirement  space  and  inequality  constraints  are  illustrated  in  Figure  2.6. 
If  a  feasible  solution  exists,  then  the  cone  generated  by  ai,a2,...,a„  must  intersect 
the  collection  of  vectors  that  are  less  than  or  equal  to  the  requirement  vector  b. 
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Figure  2.6  Requirement  space  and  inequality  constraints:  (a)  System  1  is  feasible, 
(b)  System  2  is  infeasible. 


2.1.S  Optimality.  Given  a  basic  feasible  solution  xb  =  B“^b,  where  B  is 
a  nonsingular  m  x  m  submatrix  of  A,  with  zq  —  cbXb  for  the  LP  Ax  =  b,  x  >  0, 
min  z  =  cx.  If  Zj  —  Cj  <  0  for  every  column  in  A  where  Zj  =  then  zq  is 

the  minimum  value  of  subject  to  the  constraints,  and  the  basic  feasible  solution  is 
an  optimal  basic  feasible  solution  (17:97). 

2.1.4  Critical  Regions.  The  set  of  RHS  vectors  satisfying  the  feasibility 
condition  is  called  the  critical  region  for  B  where  B  is  an  optimal  basis  for  all  RHS 
vectors  b  such  that  B“^b  >  0.  Each  critical  region  is  a  convex  cone,  the  number  of 
critical  regions  associated  with  an  LP  problem  is  finite,  and  POS(A)  is  the  union  of 
the  critical  regions  (21:11). 

2.1.5  Objective  Function  Value  as  a  Function  of  the  RHS  Vector.  Let  $(b) 
denote  the  optimal  objective  function  value  as  a  function  of  the  RHS  vector  b.  For 
each  b,  either  1)  there  exists  an  optimal  basis  B*  for  every  critical  region  z  =  1, ....,  n 
such  that  $(b)  =  CbXb  =  CB[B*]“^b  or  2)  b  ^  POS{A)  and  the  feasible  region 
is  empty  (21:15).  $  is  piecewise  linear  over  POS(A).  Since  #  is  a  continuous  and 
concave  function  over  POS(A)  (see  Murty,  Theorems  8.4  and  8.5,  p.  289)  (29),  by 
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the  Weierstrauss  Approximation  Theorem,  ^  can  be  approximated  by  a  polynomial 
(2:481). 

2.2  Response  Surface  Methodology 

Response  surface  methodology  (RSM)  was  introduced  by  G.E.P.  Box  and  K. 
B.  Wilson  in  1951  (7).  It  was  developed  for  the  purpose  of  approximating  a  surface 
in  a  region  of  interest  in  order  to  search  for  an  optimal  value. 

RSM  comprises  a  set  of  statistical  and  mathematical  techniques  for  empirical 
model  building  and  exploitation  that  encompasses  1)  Designing  a  series  of  exper¬ 
iments  that  will  yield  adequate  and  reliable  measurements  of  the  response(s);  2) 
Analyzing  the  results  of  those  experiments  to  determine  a  mathematical  model  that 
best  fits  the  data  collected,  and  3)  Searching  for  the  optimal  settings  of  the  input 
variables  that  produce  a  desired  response  (25). 

To  locate  the  combination  of  factors  which  provide  the  extreme  value  on  an 
implicit  or  explicit  response  surface,  direct  methods  and  approximation  methods 
are  used.  The  direct  methods  consist  of  direct  elimination  techniques,  hill  climbing 
techniques,  and  the  single-factor  method.  In  contrast  to  direct  methods  which  use  a 
sequential  method,  approximating  methods  consist  of  a  single  large  experiment.  The 
data  points  are  chosen  randomly  or  by  the  factorial  approach  of  selecting  all  possible 
combinations  of  the  factor  levels.  Then,  the  response  is  evaluated  at  each  of  these 
points  in  order  to  approximate  the  true  functional  relationship  between  the  response 
and  the  factors.  This  approximation  is  usually  developed  by  fitting  a  polynomial  to 
the  data.  Then,  this  polynomial  is  manipulated  by  classical  methods  to  approximate 
the  true  optimal  response  (31:29). 

2.2.1  Experimental  Design.  Box  (5:31)  defines  an  experimental  design  as 
“Each  combination  of  levels  of  the  factors  corresponding  to  a  point  in  the  factor 


2-6 


space  and  the  pattern  of  such  points  used  to  elucidate  the  [response]  surface  is  called 
the  experimental  design.” 

Montgomery  (27:27)  defines  experimental  design  as  “a  process  of  inducing  pur¬ 
poseful  changes  in  the  input  variable(s)  in  order  to  observe  and  model  the  changes 
in  the  response(s).” 

The  general  problem  of  experimental  design  in  response  surface  methodology  is 
to  choose  a  design  such  that  the  assumed  polynomial  that  is  fitted  by  the  least  squares 
regression  most  closely  represents  the  true  function  over  the  factor  space  of  interest. 
Experimental  designs  are  selected  so  that  the  experimental  error  is  minimized  or  the 
bias  due  to  the  model  misspecification  is  minimized. 

Today,  a  class  of  experimental  designs  known  as  factorial  designs  are  widely 
used.  Factorial  designs;  1)  allow  many  simultaneous  comparisons  to  be  made;  2) 
yield  highly  efficient  parameter  estimates;  and  3)  are  computationally  simple  to 
analyze  (6:106).  A  factor  is  a  variable  or  characteristic  that  is  changed  during  the 
experimentation.  Factors  may  have  different  levels  which  are  the  values  to  be  set 
during  the  course  of  experimentation.  A  specified  combination  of  levels  for  all  factors 
included  is  called  a  design  point.  So,  an  experimental  design  can  also  be  defined  as 
a  schedule  of  the  design  points  to  be  investigated  in  an  experiment. 

Experimental  designs  are  usually  divided  into  classes  according  to  the  degree 
(order)  of  polynomial  to  be  produced.  Examples  are  first  order  designs  and  second 
order  designs.  Furthermore,  factorial  designs  can  be  grouped  according  to  the  design 
points  considered.  A  complete  (or  full)  factorial  design  is  a  design  containing  every 
possible  combination  of  the  levels  of  all  factors.  If  there  are  k  factors  with  each  factor 
i  =  1,2,  ...,fc  having  n;  levels,  then  the  resulting  full  factorial  design  contains 
rai  X  ^2  X  ...  X  Uk  design  points;  if  ni  =  n2  =  ...  =  Uk  =  n,  then  it  is  referred  as  an  n* 
design.  For  example,  2“*  indicates  that  it  is  a  two  level  full  factorial  design  on  four 
factors  (variables)  and  16  design  points  (trials). 
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Fractional  factorial  designs  are  designs  which  require  only  a  fraction  of  the 
number  of  trials  required  by  a  full  factorial  design  in  order  to  provide  the  minimal 
information  requirements  (31:37).  For  example,  a  design  is  a  ^  fraction  of 
a  2^  design.  A  half  replicate  is  a  one-half  fraction  or  a  2^“^  design;  a  quarter 
replicate  is  a  one- fourth  fraction  or  a  2*^“^  design.  Although  fractional  factorial 
designs  are  very  practical  in  that  they  save  money  and  effort  for  large  experiments, 
the  main  disadvantage  of  their  use  is  that  there  is  some  confounding  of  the  effects. 
Confounding  means  that  “the  statistic  which  measures  a  main  effect  will  be  identical 
to  the  statistic  that  measures  certain  other  interaction  effects”  (19:426).  For  example, 
many  fractional  factorial  designs  are  used  to  measure  main  and  two-factor  interaction 
effects,  and  it  is  assumed  that  the  higher-order  effects  with  which  the  main  effects 
are  confounded  are  very  small. 

Resolution  is  the  length  of  the  shortest  effect  in  the  defining  relation  in  a  two 
level  fractional  design.  Let  us  denote  the  factor  by  /.  If  a  2^“^  design  is  of  resolution 
i?,  then  no  /effect  will  be  confounded  with  any  other  effect  containing  less  than  R-f 
factors.  For  example,  for  a  resolution  IV  fractional  factorial  design  of  2^"^,  no  main 
effects  are  aliased  with  any  other  main  effects  or  with  any  two-factor  interactions; 
however,  two-factor  interactions  are  aliased  with  each  other. 

2.2.2  Regression  Analysis.  Regression  analysis  was  first  developed  by  Sir 
Francis  Gallon  in  the  latter  part  of  the  19th  Century.  Gallon  had  studied  the  relation 
between  the  height  of  fathers  and  sons  and  noted  that  the  heights  of  sons  of  both 
tall  and  short  fathers  appear  to  “revert”  or  “regress”  to  the  mean  of  the  group. 
He  considered  this  tendency  to  be  a  regression  to  “mediocrity” .  Gallon  developed  a 
mathematical  description  of  this  regression  tendency  which  is  the  ancestor  of  today’s 
regression  models  (30). 

Neter,  et  al.  (30)  summarizes  the  basic  concept  as: 

A  regression  model  is  a  formal  means  of  expressing  the  two  essential 

ingredients  of  a  statistical  relation:  1)  A  tendency  of  the  dependent  vari- 
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able  Y  to  vary  with  the  independent  variable  in  a  systematic  fashion.  2) 

A  scattering  of  points  around  the  curve  of  statistical  relationship.  These 
two  characteristics  are  embodied  in  a  regression  model  by  postulating 
that:  1)  There  is  a  probability  distribution  of  Y  for  each  level  of  X.  2) 

The  means  of  these  probability  distributions  vary  in  some  systematic 
fashion  with  X. 

The  functional  relationship  between  the  responses  and  the  variables  can  be 
defined  as 


Y  =  X/3  +  e  (2.3) 

where  Y  is  s  x  1  vector  of  responses,  X  is  a  s  x  5  matrix  of  the  coded  levels  of  the 
variables,  /3  is  a  g  x  1  vector  of  parameters,  e  is  a  s  x  1  vector  of  independent  normal 
random  variables  with  expectation  E{e}  =  0,  s  is  the  number  of  data  points,  and  q 
is  the  number  of  variables. 

The  least  square  estimators  of  the  regression  model  are 

/3  =  (X'X)-^X'Y  (2.4) 

These  least  squares  estimators  are  also  “maximum  likelihood  estimators” ,  and 
they  are  “unbiased,  minimum  variance  unbiased,  consistent,  and  sufficient”  (30:238). 

The  fitted  values  and  residual  terms  are  represented  by 

Y  =  X0  (2.5) 

and 

e  =  Y -Y  =  Y -X0  (2.6) 

where  Y  is  the  s  X  1  vector  of  fitted  values,  and  e  is  a  s  X  1  vector  of  residuals. 
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The  goal  of  least  squares  regression  is  to  determine  an  estimate  /3  of  /3  so  that 
the  residual  sum  of  squares  e  e  is  minimized.  Generally,  e,  the  residuals,  are  assumed 
to  be  caused  by  random  measurement  error  (21:20). 

2.3  Kriging 

Kriging  is  a  geostatistical  estimation  technique  developed  by  Matheron  in  the 
1960’s  to  estimate  ore  reserves  from  spatial  core  samples.  It  is  named  after  D.G. 
Krige  who  was  probably  the  first  to  make  use  of  spatial  correlation  and  best  linear 
unbiased  estimator  (B.L.U.E.)  in  the  field  of  mineral  recourses  evaluation  (11:237). 

Journel  describes  kriging  as  “a  local  estimation  technique  which  provides  the 
best  linear  unbiased  estimator  of  the  unknown  characteristics  studied”  (23:304). 
So,  kriging  is  associated  with  “B.L.U.E.”.  Kriging  is  “best”  because  it  minimizes 
the  variance  of  the  errors;  it  is  “linear”  because  its  estimates  are  weighted  linear 
combinations  of  the  available  data;  it  is  “unbiased”  since  it  tries  to  have  the  mean 
residual  or  error  equal  to  zero.  The  distinguishing  feature  of  kriging  from  other 
estimation  methods  is  its  aim  of  minimizing  the  error  variance  (20:278). 

2.3.1  Kriging  Equations.  In  kriging,  the  value  of  an  unknown  point  is 
estimated  using  a  weighted  average  of  the  values  of  previously  sampled  points  in 
such  a  way  that  the  points  close  to  the  point  that  we  are  trying  to  estimate  must 
have  more  weight  then  the  ones  far  away  (11). 

The  estimator  is  given  in  matrix  form  by 


Yp  =  W'Y  (2.7) 

where  Y  is  the  estimate,  W  is  the  s  x  1  vector  of  the  weights,  and  Y  is  the  s  x  1 
vector  of  the  values  for  the  known  sample  points. 
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To  determine  the  combination  which  minimizes  the  estimation  error,  the  esti¬ 
mation  variance  must  be  defined.  The  estimation  variance  for  the  general  unbiased 
linear  estimator  is 

,7=  =  (2.8) 

!  =  1  4=1  j=l 

where  is  the  estimator  variance,  cu*  and  u>j  are  the  weights,  s  is  the  size  of  the 
sample  set,  Y  is  the  point  to  be  estimated  and  uJi^^Si^Yp)  is  known  as  the  average 
semi-variogram.  Clark  describes  a  semi-variogram  as  “a  graph  (and/or  formula) 
describing  the  expected  difference  in  value  between  pairs  of  samples  with  a  given 
relative  orientation”  (10:11). 

If  the  weights  sum  to  one  and  there  is  no  drift  (no  trend),  the  estimates  are 
considered  unbiased  (12:384).  The  weights  are  found  by  using  the  following  kriging 
equations  in  matrix  form: 

0  1  1  •••  1  A  1 

f  Tfcii  T/4i2  '  ’  ’  'fhis  'yhip 

1  7/421  7/422  •••  7/423  ^2  =  7/42P  (2-9) 

1  Ihsi  lhs2  '  '  '  'Ifhss  ‘Jhsp 

where  are  the  semi-variograms  between  the  points  in  the  neighborhood,  jhip 
are  calculated  semi-veriograms  between  the  points  in  the  neighborhood  and  the  es¬ 
timated  point,  A  is  a  Lagrange  multiplier  used  to  solve  the  simultaneous  equations, 
and  hij  are  the  distances  between  sampled  points. 

2.S.S  Types  of  Kriging.  Although  there  are  several  types  of  kriging  such 
as  lognormal,  disjunctive,  point,  block,  universal,  and  positive  kriging,  they  are  all 
related,  and  they  are  the  refined  versions  of  the  weighted  moving  average  techniques 
(18:25). 
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Point  kriging  (ordinary  kriging).  This  kriging  technique  is  used  when  we  as¬ 
sume  that  there  is  no  trend  or  drift  in  data.  To  calculate  the  point  kriging  weights, 
the  pattern  of  spatial  continuity  that  one  desires  the  random  function  model  to  ob¬ 
tain  must  be  decided.  Random  function  models  are  produced  by  variograms.  Some 
standard  models  are  generated  for  variograms.  The  smallest  standard  model  is  a 
linear  model  of  the  form  'y{h)  =  ah  +  b,  which  does  not  have  a  sill  variance  of  the 
samples.  Probably  the  most  common  standard  model  is  the  spherical  model  which 
uses  the  parameters  a,  C,  and  Cq.  The  first  parameter  a  is  called  range,  and  it  pro¬ 
vides  a  distance  beyond  which  the  variogram  or  covariance  value  remains  essentially 
constant.  Cq  is  called  the  nugget  effect,  and  it  provides  a  discontinuity  at  the  origin. 
C  is  used  in  conjunction  with  Cq  to  determine  the  sill,  (C  -|-  Cq).  The  sill  is  the 
variogram  value  for  very  large  distances,  7(00).  It  is  also  the  covariance  value  for 
I  h  I  (20:293).  The  form  of  the  spherical  model  is  given  by 


lih)  = 


c{fa-h)  +  Co 
C  +  Co 
0 


\i  h  <  a 
a  h  >  a 
ii  h  =  0 


(2.10) 


The  shape  of  this  model  is  illustrated  in  Figure  2. 7. 


Figure  2.7  Spherical  Model 
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Block  Kriging.  Block  estimate  is  an  estimate  of  the  average  value  of  a  variable 
within  the  prescribed  local  area.  So,  “one  method  for  finding  such  an  estimate  is  to 
discretize  the  local  area”  (20:323).  Since  this  method  is  simple  but  computationally 
time  consuming,  block  kriging  is  developed  as  an  alternative  to  this  method  by 
modifying  the  point  kriging  system. 

Block  kriging  uses  the  same  covariance  matrix  as  point  kriging.  However, 
the  covariance  vector  consists  of  covariance  values  between  the  random  variables  at 
the  sample  locations  that  we  are  trying  to  estimate.  For  block  estimation,  these 
covariances  are  point-to-block  covariances  while  they  are  point-to-point  covariances 
for  point  estimation.  So,  it  is  possible  to  convert  a  point  kriging  system  to  an  ordinary 
block  kriging  system  by  making  a  single  change  (20:325).  The  main  advantage  of 
block  kriging  compared  to  point  kriging  is  that  it  provides  an  estimate  of  the  block 
average  with  the  solution  of  only  one  kriging  system;  however,  block  kriging  requires 
more  computations  than  point  kriging. 

Universal  Kriging.  This  technique  is  used  when  trend  is  present.  There  are 
three  steps  in  Universal  Kriging:  1)  estimation  and  removal  of  the  drift;  2)  kriging 
of  the  nonstationary  residuals  to  obtain  needed  estimates,  and  3)  combination  of 
estimated  residuals  with  the  drift  to  obtain  estimates  of  the  actual  surface  (12:397). 

Removing  the  drift  can  be  done  by  estimating  a  polynomial  for  the  drift  and 
then  subtracting  this  drift  from  the  data.  However,  combining  the  neighborhood  size 
with  the  drift  equations  makes  this  technique  difficult  to  apply  (16:259). 

2.3.3  Conclusion. 

This  chapter  provides  a  general  review  of  the  literature  of  linear  programming 
and  the  basic  aspects  of  techniques  that  are  used  by  Johnson,  et  al.  to  create  a 
methodology  for  performing  optimality  analysis.  For  optimality  analysis,  metamod¬ 
els  are  developed  by  using  experimental  design  and  least  squares  regression.  By 
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applying  one  of  the  kriging  methods,  these  metamodels  can  be  improved  to  get  more 
accurate  estimates. 
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III.  Methodology 


If  the  optimal  objective  function  value  is  defined  as  the  response,  and  the  levels 
of  RHS  vectors  are  defined  as  predictors,  RSM  can  be  applied  to  predict  the  optimal 
objective  function  value  based  on  the  values  of  the  elements  of  the  RHS  vector. 

The  response  surface  methodology  uses  experimental  design  and  least  squares 
regression  to  develop  a  simple  polynomial  model.  Since  this  model  is  produced 
from  another  model,  in  our  case  from  a  large  LP,  it  is  called  a  metamodel  (21:14). 
Metamodels  are  used  to  predict  the  optimal  objective  function  value  based  on  the 
values  of  the  elements  of  the  RHS  vector  since  the  optimal  objective  function  value 
is  defined  as  a  function  of  the  RHS  vectors  in  the  metamodels. 

Figure  3.1  shows  the  metamodels  for  STORM  that  were  developed  in  this 
thesis.  The  first  metamodel  describes  the  relationship  between  the  objective  function 
value,  total  flying  cost,  and  the  RHS  vectors,  tonnage  requirement  (TREQ)  and 
frequency  requirement  (FREQ).  The  first  factor,  TREQ,  is  a  382  x  1  vector  and  the 
elements  of  the  vector  vary  from  0.01  to  857  in  tons.  The  second  factor,  FREQ,  is 
a  151  X  1  vector  whose  elements  range  from  4.29  to  24.5  in  visits  per  month  between 
city  pairs. 


Table  3.1  Four  Areas  in  the  AMC  Channel  System 


Area 

Region 

1 

Central  America,  Hawaii,  Guam,  Caribbean,  Alaska 

2 

Germany,  England,  Azores 

3 

Okinawa,  Japan,  South  Korea,  Italy,  Spain 

4 

North  Atlantic,  Singapore,  Australia,  New  Zealand, 
Saudi  Arabia,  Turkey,  Israel 
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Metamodels  2,  3,  and  4  show  the  relationships  between  the  total  cost  and  the 
tonnage  and  frequency  requirements  for  four  areas  designated  by  AMC.  Table  3.1 
shows  these  areas.  The  Al,  A2,  A3,  and  A4  vectors  were  simply  derived  from 
TREQ  and  FREQ  RHS  vectors  via  separating  the  elements  of  TREQ  and  FREQ 
vectors  into  the  related  areas.  Thus,  the  sum  of  the  number  of  elements  in  each  area 
vector  is  equal  to  the  number  of  elements  in  the  TREQ  and  FREQ  vectors. 

Finally,  metamodels  5  and  6  describe  the  relationship  between  the  total  flying 
cost  and  the  region  of  interest  in  the  dominating  area  associated  with  TREQ  and 
FREQ.  Region  vectors  were  obtained  from  the  area  vector  of  interest  in  the  same 
way  as  the  area  vectors  obtained  from  the  TREQ  and  FREQ  vectors. 


Figure  3.1  Research  Directions 
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The  following  sections  discuss  how  the  techniques  used  in  the  methodology  of 
Johnson,  et  al.  as  illustrated  in  Figure  3.2  are  applied. 


A 

Y=ZP 


Figure  3.2  The  Methodology  of  Johnson, ei  al.  for  Linear  Programs 
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3.1  Experimental  Design  for  Linear  Programming 

As  stated  in  the  second  chapter,  experimental  designs  are  selected  so  that  the 
experimental  error  variance  is  minimized.  In  our  study,  we  obtain  our  response 
by  solving  STORM  for  each  design  point.  In  this  case,  no  experimental  error  is 
associated  with  the  solution  from  the  LP  model.  Since  there  is  no  experimental 
error,  the  purpose  is  to  minimize  the  misspecification  bias  in  this  research. 

Minimization  of  the  bias  using  minimum  bias  estimation  was  first  suggested  by 
Karson,  Manson,  and  Hader  (24).  Karson,  et  al.  uses  an  estimator  which  achieves 
minimum  integrated  squared  bias.  This  estimator  does  not  depend  on  the  unknown 
coefficients,  and  it  attains  the  same  minimum  integrated  squared  bias  whatever  the 
design  is.  When  the  ’’all  bias”  designs  are  used  in  experimental  design,  the  minimum 
bias  estimation  of  Karson,  et  al.  is  equal  to  the  least  squares  estimation  (6). 

In  this  research,  two  experimental  designs  are  used.  The  first  design  develops 
a  data  set  containing  the  RHS  vectors  for  the  LP  and  the  corresponding  optimal 
objective  function  values  obtained  by  solving  the  LP  for  each  design  point.  This 
data  set  is  called  the  test  set.  Metamodels  are  developed  by  applying  least  squares 
regression  to  the  test  set.  The  second  experimental  design  is  used  to  obtain  a  data 
set  which  is  called  the  validation  set.  The  validation  set  is  chosen  from  the  interior 
of  the  test  set.  Metamodels  are  validated  by  the  second  experimental  design  via 
comparing  the  optimal  objective  function  values  associated  with  the  validation  set  to 
the  predictions  for  the  optimal  objective  function  values  obtained  from  metamodels. 

The  first  step  in  the  experimental  design  is  to  identify  the  admissible  set  of 
RHS  vectors.  Since  the  metamodels  will  only  be  valid  over  this  admissible  set  of 
RHS  vectors,  it  must  be  defined  such  that  all  RHS  vectors  of  potential  interest  are 
included.  In  this  study,  admissible  sets  are  constructed  by  increasing  or  decreasing 
the  RHS  vectors  of  interest  by  a  prespecified  percentage  (usually  10  percent). 
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Table  3..2  Experimental  Designs  Used  to  Develop  Metamodels  In  This  Research 


Metamodel 

Number  of  Factors 

Factorial  Design 

1 

2  TREQ,  FREQ 

2'^  full 

2 

4  TREQ  (A1,A2,A3,A4) 

fractional 

3 

4  FREQ  (A1,A2,A3,A4) 

4 

8  TREQ(A1-A4),FREQ(A1-A4) 

5 

5  TREQ(Al:CAm,Ha,Gu,Ca,Al) 

6 

6  FREQ(A3:Ja,Ks,It,Sp) 

fractional 

The  next  step  in  the  experimental  design  is  to  code  the  level  of  the  constraints. 
Coding  is  achieved  using  a  simple  transformation  given  by 


(3.1) 


where  bi  is  the  actual  numerical  level  of  the  ith  constraint.  In  this  study,  since  we 
use  two-level  factorial  designs,  we  can  say  that  bi  has  bi^max  and  bi,min  indicating 
the  higher  and  the  lower  levels  of  the  constraints.  Here  bi_o  is  fhe  midpoint  between 
bijinax  and  bi^min.  If  Sj  is  taken  to  be  (bj^max  bi^inin)/2,  then  bi^max  and  bj^min  are 
mapped  to  1  and  -1,  respectively. 

The  experimental  designs  used  to  develop  the  metamodels  of  interest  in  this 
research  were  determined  according  to  the  number  of  RHS  vectors  they  were  dealing 
with.  These  designs  were  derived  from  the  tables  by  Box  and  Draper  (5:164-165). 

For  example,  to  create  a  model  which  represents  the  relationship  between  four 
areas  and  the  total  flying  cost  due  to  tonnage  requirement  or  frequency  requirement, 
a  2|y^  design  was  used  since  there  are  four  factors  associated  with  TREQ  or  FREQ. 
The  designs  used  in  this  research  are  shown  in  Table  3.2.  The  construction  of  these 
designs  is  described  in  the  next  chapter. 
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3.1.1  Least  Squares  Regression.  In  this  research,  we  approximate  the  LP 
model,  specifically  STORM,  with  a  simple  model.  To  achieve  this,  we  apply  least 
squares  regression  to  the  data  obtained  from  the  experimental  design  phase  of  the 
methodology  of  Johnson,  et  a/.Least  squares  regression  develops  an  initial  metamodel 
approximating  $  with  a  first  order  or  second  order  polynomial  since  those  types  of 
polynomials  are  able  to  define  some  concave  surfaces  such  as  hyperplanes  and  dome¬ 
shaped  surfaces  (21:19). 

Generally,  we  wish  to  elucidate  a  model 

y  =  /(«,«)  +  t  (3-2) 

where 

E(y)  =  y  =  /(«,«)  (3.3) 

is  the  mean  level  of  the  response  y  which  is  affected  by  k  variables  (^i,  ^2?  =^*- 

The  model  also  consists  of  q  parameters  {6i,02,  ..,0g)  =6'  and  e  which  is  an  exper¬ 
imental  error  due  to  the  misspecification.  To  examine  this  model,  first  we  make  a 
series  of  experimental  runs  with  s  different  sets  of  conditions,  ^1,^2?  ^s-  Thus,  we 

choose  the  design  points,  and  then,  we  solve  the  LP  model  for  each  design  point  to 
get  the  response  values  of  yi,y2,  ..,ys- 

If  we  compute  for  each  of  the  experimental  runs  ^1,^2?  then  we 

obtain  s  discrepancies  {2/1 {Vs- Least  squares 
regression  selects  the  value  that  makes  the  sum  of  squares  of  these  discrepancies  as 
small  as  possible.  Thus,  it  tries  to  find  the  best  estimate  of  0  in  the  form  of 

S{e)  =  t,{y.  -  (3.4) 

where  S(0)  is  called  the  sum  of  squares  function.  Note  that,  for  any  specific  selection 
of  one  of  the  q  parameters  in  6,  there  is  a  specific  value  of  S(0).  The  minimizing 
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selection  of  9  is  called  the  least  squares  estimate,  and  it  is  denoted  by  0.  In  general, 
the  goodness  of  the  least  squares  estimates  of  6  depends  on  the  distribution  of  the 
errors  (6:35).  For  that  reason,  as  stated  in  Chapter  2,  the  least  squares  estimates  are 
appropriate  if  it  is  assumed  that  the  experimental  errors  =  Vu  —  yu,u  =  1,2,  ..,s 
are  statistically  independent,  with  constant  variance,  and  normally  distributed. 

The  computation  of  least  squares  estimates  is  very  simple  when  the  response 
function  is  linear  in  the  parameters,  with  the  form 


y  —  f{it9)  —  9\Z\  +  02-^2  +  9qZq 


(3.5) 


where  the  zs  are  known  constants.  In  the  experimental  design,  zs  are  known  functions 
of  the  experimental  conditions(runs),  ^25  Equivalently  we  can  say  that  they 
are  the  coded  forms  Xi  =  {(i  -  (i,o)/Si,i  =  1,2,  where  (ifi  and  5,  are  suitable 
location  and  scale  factors,  respectively.  Adding  the  experimental  error  e  =  y  —  y,vre 
have  the  model 


y  —  +  ^2^2  + . +  9qZq  +  e  (3-6) 


Box  and  Draper  call  this  the  linear  model  since  it  is  linear  in  the  parameters. 

If  Equation  (3.5)  has  been  postulated  and  experimental  conditions  {ziu,  Z2u,  Zqu), 

u  =  1, 2,  ..,s  have  been  run,  yielding  ?/i, y2,  ?/«,  have  been  obtained,  then  the  model 

relates  the  observations  to  the  known  ZiJs  and  the  unknown  0,’s  by  s  equations 
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yi  =  Oizu  +  02^21  +  — 

...  -f  0qZq\  +  Ci 

(3.7) 

y2  —  0lZi2  +  02^22  +  ••• 

...  +  0qZq2  +  e.2 

(3.8) 

Vs  —  01^13  +  02^23  +  ••• 

...  -|-  0qZqs  +  eg 

(3.9) 

These  equations  can  be  written  in  matrix  form  as 

y  =  Z0  +  e  (3.10) 

where  Y  is  an  5  x  1  vector  of  responses,  Z  is  an  s  X  ^  matrix  of  known  constants,  6 
is  the  q  X  1  vector  of  unknown  coefficients,  and  e  is  an  s  X  1  vector  of  experimental 
error. 

As  we  can  see,  the  experimental  design  determines  the  Z  matrix.  To  estimate 
the  q  parameters  of  the  model  separately,  we  should  employ  an  experimental  design 
that  provides  a  s  x  q  matrix  with  full  column  rank  with  Z’Z  nonsingular,  so  the 
solution  to  equations  (3.7)  through  (3.9)  can  be  written 

e  =  (Z'Z)-iZ'y,  (3.11) 


and  it  is  called  the  least  squares  estimator. 


3.1.2  Kriging  of  the  Residuals.  It  is  possible  to  improve  the  metamodel 
obtained  from  the  least  squares  regression  by  kriging.  In  this  study,  simple  point 
kriging  is  applied  to  the  residuals  that  are  obtained  from  the  test  metamodel.  Thus, 
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in  Equation  (2.7),  Y  is  replaced  by  e,  and  Yp  is  replaced  by  Cp  where  e  is  a  5  x  1 
vector  of  estimated  residuals  from  the  regression  metamodel,  and  ip  is  the  estimate 
of  the  regression  residual  to  be  obtained  by  kriging.  So,  the  equation  for  estimation 
of  the  error  can  be  written  as 


4  =  W'e  (3.12) 

The  distance  matrix  is  constructed  by  calculating  the  Euclidean  distance  be¬ 
tween  each  design  point.  The  covariance  structure  of  the  points  are  derived  by  using 
one  of  the  variogram  models  summarized  in  Section  3.2.  For  example,  the  spherical 
model  can  be  applied  to  the  distance  matrix.  By  treating  the  residuals  from  the 
test  metamodel  as  known  points,  the  kriged  residuals  are  estimated  by  finding  the 
weights  from  Equation  (2.9),  and  then  applying  these  weights  to  Equation  (3.13). 

Note  that  the  performance  of  point  kriging  depends  on  the  covariance  structure 
of  the  residuals.  The  values  of  the  parameters  that  are  used  in  the  selected  variograms 
i.e.  Co,  C,  and  a  in  the  spherical  model,  must  be  chosen  so  that  the  error  is  minimum. 

3.1.3  Combination  of  Regression  and  Kriging.  The  estimated  residuals 
obtained  by  kriging  are  simply  added  to  the  regression  metamodels  to  get  more  ac¬ 
curate  metamodels.  Thus,  we  decrease  the  error  between  the  true  objective  function 
value  obtained  by  solving  the  LP  and  the  objective  function  value  predicted  by  the 
metamodel. 

3.1.4  Measuring  the  Performance  of  the  Metamodels.  The  predictions 
obtained  from  the  regression  metamodel  and  from  the  metamodel  developed  by  using 
kriging  in  conjunction  with  regression  are  compared  to  the  true  optimal  objective 
function  values  of  the  validation  design  by  using  the  measures  MSE,  MAPE,  and 
PC. 


3-9 


Mean  squared  error  (MSE)  is  a  measure  of  accuracy  computed  by  squaring  the 
individual  error  for  each  item  in  the  data  set,  and  then,  finding  the  average  or  mean 
value  of  the  sum  of  those  squares.  The  mean  squared  error  gives  greater  weight  to 
large  errors  than  to  small  errors  since  the  errors  are  squared  before  being  summed. 
MSE  is  defined  by 


MSE  = 


EUAYi  -  Yif 


(3.13) 


and  it  is  an  unbiased  estimator  of  for  the  regression  model. 

Mean  absolute  percentage  error  (MAPE)  is  the  mean  or  average  of  the  sum  of 
all  of  the  percentage  errors  for  a  given  data  set  taken  without  regard  to  sign.  Thus, 
their  absolute  values  are  summed  and  the  average  computed.  MAPE  is  given  as 
follows; 


MAPE  =  -  2 


^r=i 


Yi-Yi 


X  100, 


(3.14) 


Percentage  of  predictions  closer  to  the  true  optimal  solution  (PC)  is  defined 


by 


PC  =  -  X  Number  of  Predictions  Closer  x  100  (3.15) 

s 

PC  is  used  to  compare  the  regression  metamodel  to  the  improved  metamodel  ob¬ 
tained  by  both  regression  and  kriging. 
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IV.  Results  and  Discussion 


This  chapter  presents  the  numerical  results  and  their  interpretations.  In  this 
thesis,  six  metamodels  of  STORM  are  developed  to  check  the  applicability  of  the 
methodology  of  Johnson,  et  al.  The  test  design,  validation  design,  and  kriging  results 
are  given  for  the  first  order  designs.  Then,  the  results  concerning  second  order  designs 
and  their  interpretations  follow.  The  test  metamodels  are  developed  by  perturbing 
the  RHS  vectors  by  10  percent  (A  =  ±0.1).  Refer  to  Figure  3.1  as  an  overview 
device  for  the  following  sections.  Note  that  the  metamodels  are  given  in  the  coded 
form  of  the  factors,  and  denoted  by  Zi,i  =  1,  ...,s  in  the  following  tables. 

4-1  Metamodel  1 

Metamodel  1  shows  the  effects  of  the  factors.  Tonnage  Requirement  (TREQ) 
and  Frequency  Requirement  (FREQ)  on  the  optimal  objective  function  value.  Total 
Flying  Cost  (in  dollars).  Here,  TREQ  is  a  RHS  vector  that  contains  382  elements. 
Thus,  there  are  382  city  pairs  involved  in  cargo  delivery.  FREQ  is  a  RHS  vector 
containing  151  elements.  This  means  that  151  city  pairs  have  specified  frequency 
levels  that  must  be  attained  during  a  period.  Appendix  B  presents  the  data  sets  of 
TREQ  and  FREQ. 

4.1.1  First  Order  Design  of  Metamodel  1.  Table  4O  shows  the  experi¬ 
mental  design  and  the  regression  results  of  Metamodel  1.  Notice  that  is  very 
close  to  1,  and  the  same  amount  of  residual  error  occurs  at  every  possible  condition. 
The  metamodel  is  given  in  coded  form.  We  see  that  the  coefficient  of  TREQ  is 
approximately  twice  that  of  FREQ.  Thus,  a  10%  change  in  the  tonnage  require¬ 
ment  changes  the  total  flying  cost  twice  as  much  as  a  10%  change  in  the  frequency 
requirement. 
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Table  4.1  Metamodel  1:  Test  Design  and  Regression  Results 


Table  4.2  Metamodel  1:  Validation  Design  and  Kriging  Results 


2^  Full  Factorial  Design 

run 

TREQ 

Y 

e 

V 

i  reg 

- X - 

A 

Y 

1 

-0.5 

-0.5 

36,773,800 

-26,680 

36,889,301 

36,862,621 

2 

0.5 

-0.5 

39,495,708 

26,680 

39,567,270 

39,593,950 

3 

-0.5 

0.5 

38,128,297 

26,680 

38,199,280 

38,225,960 

4 

0.5 

0.5 

40,747,355 

-26,680 

40,877,249 

40,850,569 

Table  4.3  Regression  and  Kriging  Performance  Comparison  of  Metamodel  1 


MSE 

MAPE 

PC 

Regression 

1.009x10^° 

0.2501 

50% 

Regression  +  Kriging 

0.2499 

50% 
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If  we  recall  that  the  original  objective  function  value  of  STORM  was  $38.69M, 
by  looking  at  the  regression  metamodel  in  Table  we  can  say  that  the  contribution 
of  the  tonnage  requirement  to  the  total  flying  cost  is  $2.68M  while  it  is  $1.3M  for  the 
frequency  requirement.  In  other  words,  the  tonnage  requirement  changes  the  total 
flying  cost  by  7%  while  the  frequency  requirement  changes  the  optimal  objective 
function  value  by  3.5%.  Also,  it  is  possible  to  say  that  our  regression  metamodel 
estimates  the  total  flying  cost  with  an  error  of  ±$92645  for  the  design  points  taken 
into  account. 

To  obtain  a  better  metamodel,  needless  to  say,  a  metamodel  with  less  residual 
error,  we  apply  simple  point  kriging  to  the  residuals  obtained  from  the  regression 
metamodel  (Refer  to  Appendix  C  for  the  application  of  simple  point  kriging).  Then, 
we  add  these  kriged  residual  estimates  (e)  to  the  estimated  optimal  objective  function 
values  Y  of  the  validation  design  that  are  estimated  by  substituting  the  validation 
design  points  into  the  regression  metamodel.  The  resulting  objective  function  value 
vector  is  denoted  by  Y  in  Table  4-2-  Metamodel  1  was  validated  by  a  2^  full  factorial 
design  where  the  factors  were  coded  at  the  (±0.5,  -0.5)  level  (Perturbed  by  5%). 
Figure  4-i  shows  the  test  and  validation  design  points  on  a  coordinate  axis.  Note 
that  the  point  (0,0)  is  the  center  point  for  both  designs. 
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Figure  4.1  Test  and  Validation  Designs  of  Metamodel  1 
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The  comparison  of  the  regression  metamodel  and  the  metamodel  obtained  by 
regression  and  kriging  is  given  in  Table  4-3.  When  we  compare  Y  and  Y  to  the  true 
optimal  objective  function  value  vector,  Y,  we  see  that  Y  of  the  first  and  the  fourth 
design  points  are  closer  to  Y,  while  Y  of  the  second  and  the  third  design  points  are 
closer  to  Y.  Thus,  the  PC  is  50%  for  both  metamodels.  Note  that  MAPE  of  the 
metamodels  are  very  close  to  each  other  indicating  that  simple  point  kriging  does 
not  improve  our  regression  metamodel  much.  The  response  surface  is  flat  in  the 
regression  metamodel,  and  the  metamodel  explains  our  system  with  a  small  error 
so  that  we  do  not  need  to  improve  the  metamodel  furthermore.  For  that  reason, 
kriging  was  not  applied  to  the  following  metamodels. 

4.1.2  Second  Order  Design  of  Metamodel  1.  The  second  order 
design  of  Metamodel  1  is  given  in  Table  4-4-  The  second  order  design  with  an 
interaction  term  fits  the  data  without  any  type  of  error  with  one  degree  of  freedom. 
The  negative  interaction  term  indicates  that  if  we  use  the  channelling  cargo  system 
effectively,  it  reduces  the  total  flying  cost  by  $92645  when  the  product  of  the  tonnage 
and  frequency  requirements  is  one.  When  we  examine  STORM,  we  see  that  attaining 
frequency  requirements  without  any  cargo  delivery  causes  an  additional  increase  in 
the  total  flying  cost.  Thus,  an  effective  channelling  system  rejects  empty  aircraft 
that  have  to  be  used  to  achieve  the  minimum  required  frequency  level  for  a  given 
period. 


Table  4.4  Second  Order  Design  of  Metamodel  1 


Y  =  38,786,290  +  2,670,483zi  +  1,303,072z2  -92,645ziZ2 

Regression 

df 

Type  I  Sum  of  Squares 

B? 

Linear 

2 

3.556x10^=" 

0.9990 

Quadratic 

0 

0 

0.0000 

Crossproduct 

1 

3.433  xlO^" 

0.0010 

Total  Regression 

3 

3.558x10^^ 

1.0000 
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4-2  Metamodel  2 


Metamodel  2  investigates  the  effects  of  the  tonnage  requirement  of  four  areas 
on  the  total  cost.  These  areas  are  abbreviated  as  Al,  A2,  A3,  and  A4,  respectively, 
where  Al  is  a  94  x  1  RHS  vector,  A2  is  a  39  x  1  RHS  vector,  A3  is  a  119  x  1  RHS 
vector,  and  finally,  A4  is  a  130  x  1  RHS  vector.  Note  that  TREQ  is  composed  of 
the  four  area  vectors. 

4.2.1  First  Order  Design  of  Metamodel  2.  When  we  look  at  the  Table 
4.5,  we  observe  that  B4  is  almost  1  showing  that  we  have  a  good  first  order  fit  to 

A 

the  data.  Also,  the  predicted  optimal  objective  function  values  Y  are  close  to  the 
true  optimal  objective  function  values  Y.  The  regression  metamodel  tells  us  that  in 
case  of  a  possible  10%  change  in  tonnage  requirement.  Area  1  and  Area  4  are  the 
dominating  factors  for  the  total  flying  cost.  Especially,  Al  has  the  greatest  effect  on 
the  total  flying  cost  by  2.12%  for  the  tonnage  requirement.  But,  we  should  keep  the 
dimensions  of  the  area  vectors  in  mind.  For  example,  although  A2  seems  the  least 
important  factor  of  all,  it  consists  of  only  39  city  pairs.  So,  we  can  say  that  the  those 
routes  actually  have  great  importance  in  cargo  delivery. 

4.2.2  Second  Order  Design  of  Metamodel  2.  Table  4-6  summarizes  the 
findings.  We  have  a  perfect  fit  of  data  as  in  Metamodel  1.  Notice  that  Area  1  has  an 
interaction  with  the  other  areas  although  the  interaction  effects  are  small.  If  we  look 
at  the  map  given  in  Appendix  D  which  shows  the  channel  cargo  system  of  AMC,  we 
see  that  Area  1  is  generally  the  first  leg  for  cargo  movement  to  the  the  other  areas. 
This  causes  the  interaction  of  Area  1  with  the  other  areas. 
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Table  4.5  Metamodel  2:  First  Order  Design  and  Regression  Results 


IV  Fractional  Factorial  Design  of  TREQ 


A2  A3  A4 


-1  -1  -1  36,159,440 

-1  -1  1  39,495,836 

1  -1  1  38,748,515 

1  -1  -1  38,853,288 

-1  1  1  38,945,744 

-1  1  -1  38,996,404 

1  1  -1  38,290,836 

1  1  1  41,525,760 


36,193,818 

39.479.478 
38,753,033 
38,830,750 
38,923,206 
39,000,922 

38.274.478 
41,560,138 


-34,378 

16,358.4 

-4,518.4 

22,538.1 

22,538.1 

-4,518.4 

16,358.4 

-34,378 


Source 

df 

Sum  of  Squares 

Mean  Square 

R2 

Model 

4 

1.5159x10^^ 

3.7799x10^2 

0.9997 

Error 

3 

3,955,648,706 

1,318,549,568.7 

C  Total 


38, 876, 978  -h  840, 844zi  -b  477, 622z2  -b  562, 708z3  +  801, 986z4 


Table  4.6  Second  Order  Design  of  Metamodel  2 


Y  =  38, 870, 632  -b  840, 844zi  -b  475, 868z2  -b  578, 696z3  +  787, 752z4 

4-425. 79ziZ2  —  19,448ziZ3  —  9, 009.81ziZ4 

Regression 

df 

Type  I  Sum  of  Squares 

R2 

Linear 

4 

1.511x10^^ 

0.9998 

Quadratic 

a 

0 

Crossproduct 

11 

Total  Regression 

D 

3.5115x102^ 

1.0000 
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Table  4.7  Metamodel  3:  First  Order  Design  and  Regression  Results 


m 

4-1 

IV 

Fractional  Factorial  Design  of  FREQ 

run 

A1 

A2 

A3 

A4 

Y 

- A - 

Y 

e 

1 

-1 

-1 

-1 

-1 

37,434,501 

37,542,443 

-107,942 

2 

1 

-1 

-1 

1 

38,973,981 

38,842,009 

131,973 

3 

-1 

1 

-1 

1 

38,445,504 

38,475,241 

-29,737.5 

4 

1 

1 

-1 

-1 

38,461,651 

38,455,945 

5,706.7 

5 

-1 

-1 

1 

1 

39,241,418 

39,235,712 

5,706.7 

6 

1 

-1 

1 

-1 

39,186,678 

39,216,415 

-29,737.5 

7 

-1 

1 

1 

-1 

38,981,620 

38,849,647 

131,973 

8 

1 

1 

1 

1 

40,041,272 

40,149,213 

-107,942 

Source 

df 

Sum  of  Squares 

Mean  Square 

i?2 

Model 

4 

3.9763x10^2 

994,085,642,302 

0.9851 

Error 

3 

59,970,102,201 

19,990,034,067 

C  Total 

7 

4.0363x10^2 

Y  =  38, 845, 828  +  320, 067zi  +  136, 683z2  +  516, 919z3  +  329, 716^ 


4.3  Metamodel  8 

Metamodel  3  shows  the  key  relationships  among  the  roles  of  the  four  areas  for 
the  frequency  requirement.  Here,  Area  1  is  a  52  x  1  RHS  vector,  Area  2  is  a  12  x  1 
RHS  vector,  Area  3  is  a  41  x  1  RHS  vector,  and  Area  4  is  a  46  x  1  RHS  vector  of 
FREQ.  As  in  Metamodel  2,  the  FREQ  data  set  is  composed  of  these  areas. 

4.3.1  First  Order  Design  of  Metamodel  3.  Table  .^.7  shows  the  test  design 
and  regression  results.  Note  that  the  E?  is  again  close  to  1.0  indicating  that  the 
metamodel  fitted  the  data  very  well.  For  the  frequency  requirement,  Area  3  which 
is  composed  of  Japan,  South  Korea,  Italy,  and  Spain  has  the  greatest  effect  on  the 
total  flying  cost  by  1.3%  when  there  is  a  10%  change  in  the  frequency  requirement. 
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4-3.2  Second  Order  Design  of  Metamodel  3.  In  Table  4-3,  we  see  that 
the  second  order  metamodel  fits  the  data  perfectly  again  without  any  type  of  error 
with  three  degrees  of  freedom.  Area  1  is  the  dominating  factor  for  the  frequency 
requirement  as  in  the  first  order  design  of  Metamodel  3.  Notice  that  Area  1  has 
interactions  with  the  other  areas  like  the  second  order  metamodel  2.  However,  in 
this  case,  Area  3  plays  the  most  important  role.  Also,  the  interaction  of  Area  1  and 
Area  4  are  significant  since  the  frequency  requirement  is  satisfied  for  Australia  and 
New  Zealand  in  Area  4  via  Hawaii  which  is  in  Area  1. 


Table  4.8  Second  Order  Design  of  Metamodel  3 


Y  =  38,970,828  +  445, 067zi  4-261,  683z2  +641, 919z3  +454,  919z4 

+73,  883ziZ2  +  56, 16OZ1Z3  —  137, 015ziZ4 

Regression 

df 

Type  I  Sum  of  Squares 

B? 

Linear 

4 

7.083  xlO^'-^ 

0.97 

Quadratic 

0 

0 

0.0000 

Crossproduct 

3 

2.19x10^1 

0.03 

Total  Regression 

7 

7.302x10^2 

1.0000 

4-4  Metamodel  4 

This  metamodel  is  the  combination  of  Metamodel  1  and  Metamodel  2.  It 
yields  the  same  insights  as  Metamodels  1  and  2  since  we  get  the  same  results  from 
the  first  order  design.  They  are  shown  in  Table  4-3.  Note  that  it  is  easy  to  see  that 
Metamodels  1  and  2  were  developed  correctly. 
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Table  4.9  Metamodel  4:  First  Order  Design  and  Regression  Results 


2®  ^iv  Fractional  Factorial  Design  of  TREQ  and  FREQ  A1  A2  A3  A4 

run 

TAl  TA2  TA3  TA4  FAX  FA2  FA3  FA4 

Y 

- A 

Y 

e 

1 

-1  -1  -1  -1  -1  -1  -1  -1 

34802682 

34996611 

-193929 

2 

1-1-1-1-11  1  1 

38584729 

38616983 

-32253.6 

3 

-11-1-11-11  1 

38326536 

38199785 

126751 

4 

1  1-1-111-1-1 

38542247 

38520104 

22143.1 

5 

-1-11-1111-1 

38219106 

38125097 

94009.1 

6 

1-11-11-1-11 

39258659 

39169893 

88765.6 

7 

-11  1-1-11-11 

37986503 

38035961 

-49457.7 

8 

1  1  1-1-1-11-1 

39679532 

39735561 

-56029 

9 

-1-1-1111-11 

38177527 

38233556 

-56029 

10 

1-1-11  1-11-1 

39883699 

39933156 

-49457.7 

11 

-11-11-111-1 

38887989 

38799224 

88765.6 

12 

1  1-11-1-1-11 

39938030 

39844021 

94009.1 

13 

-1-11  1-1-111 

39471157 

39449014 

22143.1 

14 

1-11  1-11-1-1 

39896084 

39769332 

126751 

15 

-11  1  1  1  -1  -1  -1 

39319881 

39352135 

-32253.6 

16 

11111111 

42778578 

42972506 

-193929 

Source 

df 

Sum  of  Squares 

Mean  Square 

Model 

8 

3.8215x10^3 

0.9960 

Error 

7 

1550142525548 

mmmmm 

C  Total 

15 

Y  =  38, 984, 559  +  835, 636zi  +  447, 853z2  +  591, 629z3  +  809, 559z4 
+328,  720z5  +  149, 537z6  +  494, 357z7  +  330, 656z8 
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4-5  Metamodel  5 


The  following  two  metamodels  can  be  called  the  regional  search  metamodels 
since  they  are  trying  to  find  the  key  relationships  between  the  regions  and  the  total 
cost  in  the  dominating  areas,  namely  Area  1  for  the  tonnage  requirement  and  Area 
3  for  the  frequency  requirement. 

Metamodel  5  investigates  the  effect  of  Central  America,  Hawaii,  Guam,  Caribbeans, 
and  Alaska  on  the  total  flying  cost.  In  Figure  4-10i  we  see  that  Central  America 
is  able  to  change  the  objective  function  value  by  1%.  Since  there  is  a  good  fit  to 
the  data  with  =  1.0,  we  do  not  need  to  examine  the  second  order  metamodel 
although  there  is  an  insignificant  interaction  between  Hawaii  and  the  other  regions 
because  Hawaii  is  an  important  stopover  point  in  the  middle  of  the  Pacific  Ocean. 


Table  4.10  Metamodel  5:  First  Order  Design  and  Regression  Results 


2^  ^iii  Fractional  Factorial  Design  of  TREQ  Area  1 

run 

CAm  HA  GU  CAR  AL 

Y 

- * - 

Y 

e 

1 

-1  -1-1  1  1 

38,450,021 

38,450,303 

-281.3 

2 

1  -1  -1  -1  -1 

38,732,340 

38,736,397 

-4,057.2 

3 

-1  1-1-1  1 

38,360,306 

38,360,025 

281.3 

4 

1  1-11-1 

39,309,170 

39,305,112 

4,057.2 

5 

-1-11  1  -1 

38,354,151 

38,353,869 

281.3 

6 

1-11-11 

38,993,751 

38,989,694 

4,057.2 

7 

-111-1  -1 

38,263,310 

38,263,592 

-281.3 

8 

11111 

39,554,352 

39,558,409 

-4,057.2 

Source 

df 

Sum  of  Squares 

Mean  Square 

B? 

Model 

1.6547x10^2 

330937442922 

1.0000 

Error 

2 

66,160,077.761 

33,080,038.881 

C  Total 

7 

1.65475x10^2 

Y  =  38, 752, 175  +  395, 228zi  +  119,  fiOOzg  +  39, 2I6Z3  +  164, 748z4  +  87, 433z5 
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4-6  Metamodel  6 


This  metamodel  was  developed  to  observe  the  regional  effects  on  FREQ  of 
Area  3  since  this  area  is  dominant  on  the  total  flying  cost  as  we  observe  from  Meta¬ 
model  3. 

4.6.1  First  Order  Design  of  Metamodel  6.  Tables  4-il  and  ^.i^show  that 
Japan  and  Italy  are  the  dominating  regions  on  the  total  cost  by  0.65%  for  a  10% 
change  in  the  frequency  requirements.  The  map  in  Appendix  D  showing  the  AMC’s 
Cargo  Channel  System  confirms  our  statement. 


Table  4.11  Metamodel  6;  First  Order  Design  and  Regression  Results 


2^  ^  IV  Fractional  Factorial  Design  of  FREQ  Area  3 

run 

JAPAN 

SKOREA 

ITALY  SPAIN 

Y 

'  ■  '  A  '  ' 

Y 

e 

1 

-1 

-1 

-1  -1 

38,213,241 

38,229,721 

-16,479.9 

2 

1 

-1 

-1  1 

38,760,857 

38,744,381 

16,475.8 

3 

-1 

1 

-1  1 

38,284,713 

38,268,245 

16,467.7 

4 

1 

1 

-1  -1 

38,766,451 

38,782,914 

-16,463.5 

5 

-1 

-1 

1  1 

38,700,791 

38,717,254 

-16,463.5 

6 

1 

-1 

1  -1 

39,248,391 

39,231,923 

16,467.7 

7 

-1 

1 

1  -1 

38,772,263 

38,755,787 

16,475.8 

8 

1 

1 

1  1 

39,253,968 

39,270,448 

-16,479.9 

Source 

df 

Sum  of  Squares 

Mean  Square 

B? 

Model 

4 

1.0081x10^2 

252028629857 

0.9979 

Error 

3 

2,170,540,199.8 

723,513,399.93 

C  Total 

7 

1.0103x10^2 

Y  =  38, 750, 084  +  257,332zi  +  19,264z2  +  243, 769z3  +  OZ4 
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j(.6.2  Second  Order  Design  of  Metamodel  6.  In  Figure  we  see  that 
the  interaction  between  Japan  and  South  Korea  is  significant  because  these  channels 
serve  the  same  region  of  interest. 


Table  4.12  Second  Order  Design  of  Metamodel  6 


Y  =  38750084  -h  257332zi  +  19264z2  +  243769z3  +  OZ4  -  16472ziZ2 

Regression 

df 

Type  I  Sum  of  Squares 

B? 

Linear 

4 

1.008x10^'-^ 

0.9979 

Quadratic 

0 

0 

0.0000 

Crossproduct 

3 

2.17x10^^ 

0.0021 

Total  Regression 

7 

1.01x10^'-^ 

1.0000 
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V.  Conclusions  and  Recommendations 

5.1  Conclusions 

The  methodology  developed  by  Johnson,  Bauer,  Moore,  and  Grant  for  the 
optimality  analysis  of  LPs  works  very  well  for  the  large  scale  linear  programming 
models  as  was  shown  in  Chapter  4.  The  metamodels  developed  for  different  sce¬ 
narios  provided  some  useful  insights  about  the  LP  model,  namely  STORM,  and  the 
metamodels  described  the  key  relationships  between  the  optimal  objective  function 
value,  the  total  flying  cost,  and  the  right-hand-side  vectors  of  interest. 

The  simple  point  kriging  did  not  contribute  a  considerable  amount  of  improve¬ 
ment  to  the  developed  metamodels  not  because  it  did  not  work  in  general,  but 
because  of  the  structure  of  the  STORM  model.  Note  that  the  response  surface  of 
the  optimal  objective  function  value  is  a  flat  surface.  For  that  reason,  the  optimal 
objective  function  value  can  be  estimated  by  a  simple  polynomial  with  remarkable 
accuracy  as  we  have  seen  in  our  study. 

Although  the  metamodels  we  have  developed  may  suggest  a  change  in  the 
optimal  routing  system  to  minimize  the  total  flying  cost  when  there  is  a  10%  change 
in  the  tonnage  and  frequency  requirements,  we  saw  that  the  gain  from  changing  the 
resources  by  10%  or  from  making  some  changes  in  the  routing  system  is  insignificant 
if  we  compare  the  results  with  the  total  flying  cost  which  is  almost  $40M. 

5.2  Recommendations 

5.2. 1  Recommendations  for  AMC.  This  thesis  effort  limited  itself  by  chang¬ 
ing  the  RHS  vectors  of  STORM  by  10%.  One  may  study  perturbing  the  RHS  vectors 
of  interest  for  various  levels  of  change.  To  give  an  incentive,  the  effect  of  the  frequency 
requirement  (FREQ)  on  the  total  flying  cost  is  summarized  in  Table  5.1  when  the 
change  in  TREQ  is  kept  at  constant  10%  level  while  the  FREQ  is  changed  by  50%, 
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Table  5.1  Effect  of  TREQ  on  the  Total  Cost  for  Various  Perturbations 


Effect  on  the  Optimal  Objective  Function  Value 

mm 

In  Million  Dollars 

In  Percentage 

■giMll 

6.47 

16.7% 

10.72 

27.7% 

±1.00 

14.36 

37.1% 

80%,  and  100%  (100%  change  means  the  values  of  TREQ  are  doubled  or  set  to 
zero) . 


5.2.2  Recommendations  About  the  Application  of  the  Methodology. 

1.  Excel  5.0  was  used  to  change  the  levels  of  the  RHS  data  set.  Although  it  was 
very  efficient,  one  may  use  the  Macro  utilities  of  the  related  software  to  make 
the  data  manipulations  in  a  shorter  time. 

2.  Simple  point  kriging  was  applied  to  the  regression  metamodel  residuals  in  this 
study.  Other  types  of  kriging  techniques  presented  in  Chapter  2,  for  example 
universal  kriging,  can  be  used  in  further  studies. 

3.  During  the  application  of  kriging,  other  types  of  models  can  be  used  rather 
than  the  spherical  model  to  get  better  estimates. 

5.3  Further  Study 

The  methodology  of  Johnson,  et  al.  can  also  be  applied  to  the  objective  function 
coefficients  of  LPs  since  the  objective  function  values  are  simply  the  duals  of  the  RHS 
vector  values  in  a  linear  programming  model. 
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Appendix  A.  STORM 

This  information  about  STORM  was  provided  by  Mr.  Alan  Whisman  of  the 
Command  Analysis  Group  at  AMC. 

The  Data  Requirements  for  the  STORM  Model 


1.  Aircraft  Routes  and  Flying  Times:  Possible  cargo  routes  must  be  specified  to 
the  STORM  model  in  advance,  along  with  flying  times  for  each  type  of  plane 
on  each  route. 

2.  Aircraft  Data:  Information  specific  to  each  aircraft  type,  such  as  payload, 
total  number  of  flying  hours  available,  and  cost  per  operating  hour.  Aircraft 
include  commercial  cargo  planes  (Boeing  747,  DC-8,  and  DC-10)  as  well  as 
AMC  military  aircraft  (C-5,  C-141,  and  C-130). 

3.  Cargo  and  Frequency  Requirements:  The  number  of  missions  to  be  flown 
is  based  upon  the  need  to  move  aircraft  and  cargo  between  certain  pairs  of 
bases.  Cargo  tonnage,  frequency  requirements,  and  possible  transshipment 
points  must  be  specified  for  each  pair  of  bases. 

4.  Cargo  Handling  Costs:  The  STORM  model  allows  cargo  to  move  either  by 
a  direct  flight  or  utilizing  one  transshipment.  Each  time  cargo  is  loaded  or 
unloaded,  a  handling  cost  of  $176  per  ton  is  incurred. 

5.  Base  Operating  Limits:  If  any  base  in  the  system  has  restrictions  on  the  max¬ 
imum  number  of  channel  flights  it  can  handle  in  a  month,  this  value  must  be 
specified. 

The  STORM  model  uses  the  preceding  information  to  solve  a  linear  program¬ 
ming  problem  which  can  best  be  described  as  follows: 

MINIMIZE:  Total  costs  incurred  for  aircraft  operating  hours  and  cargo  han¬ 
dling. 
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SUBJECT  TO: 


1.  Meeting  as  many  of  the  cargo  movement  requirements  as  possible. 

2.  Meeting  all  requirements  for  service  frequency. 

3.  Operating  within  each  aircraft  type’s  flying  hour  and  payload  limits. 

4.  Operating  within  each  base’s  limit  on  monthly  sorties. 

Outputs  of  the  model  include; 

1.  The  number  of  missions  flown  on  each  route  by  each  aircraft  type. 

2.  The  number  of  tons  that  could  not  be  delivered  for  each  origin/destination 
pair.  Since  STORM  allows  only  one  transshipment  for  each  piece  of  cargo, 
many  of  these  undelivered  requirements  can  be  met  with  two  transshipments. 

3.  For  each  movement  requirement,  the  tonnage  delivered  directly  by  a  single 
mission. 

4.  For  each  movement  requirement,  the  tonnage  transshipped  through  each  pos¬ 
sible  transshipment  point. 

5.  The  percentage  of  available  cargo  space  being  utilized  on  each  leg  of  each  route. 

The  linear  programming  formulation  used  in  STORM  has  been  implemented 
in  the  LP  modeling  language  GAMS.  GAMS  is  available  for  a  wide  variety  of  LP 
solvers  and  machines,  and  greatly  decreases  the  time  involved  in  formulating,  testing, 
and  debugging  models,  as  well  as  providing  an  excellent  report  writing  capability. 
The  description  given  below  is  designed  to  aid  in  interpreting  the  GAMS  program 
listing. 

There  are  six  types  of  variables  that  are  used  in  the  STORM  LP  model  of  the 
channel  cargo  system.  They  are  as  follow: 

1.  Y(c,cd):  The  number  of  tons  of  cargo  from  source  base  ‘c’  which  cannot  be 
delivered  to  final  destination  bases  ‘cd’. 
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2.  X(routes,  aircraft):  The  number  of  missions  to  be  flown  on  a  route  by  aircraft 
type. 

3.  D(routes,  stp,  st2):  The  number  of  tons  of  cargo  to  be  delivered  directly  from 
stop  number  ‘stp’  of  the  route  (its  origin)  to  stop  number  ‘st2’  of  the  route  (its 
destination). 

4.  S(routes,  stp,  st2):  The  number  of  tons  picked  up  at  a  transshipment  point 
(stop  number  ‘stp’  of  the  route)  and  dropped  at  its  final  destination  (stop 
number  ‘st2’). 

5.  T(routes,  stp,  st2,  cd);  The  number  of  tons  picked  up  at  its  origin  (stop  ‘stp’) 
and  dropped  off  at  a  transshipment  point  (stop  ‘st2’)  for  eventual  delivery  to 
its  final  destination  base  ‘cd’. 

6.  SLK(routes,  stp):  The  unused  cargo  capacity  on  the  ‘stp’th  leg  of  the  route. 

Generating  variables  for  all  combinations  of  routes,  bases  and  aircraft  types 
would  yield  an  unmanageable  number  of  redundant  and  unnecessary  variables  in  the 
LP  formulation.  To  prevent  this,  a  FORTRAN  program  examines  the  set  of  routes, 
transshipment  points,  and  cargo  requirements  and  builds  only  those  combinations 
required  for  the  current  data.  The  FORTRAN  program  builds  a  number  of  data  sets 
which  are  used  as  “INCLUDE”  input  files  in  the  GAMS  program.  The  created  files 
have  extensions  .TMP  and  .VAR. 

There  are  nine  types  of  constraints  in  the  model.  They  are  listed  below  by 
their  GAMS  names  and  described  briefly: 

1.  TREQ(c,  cd):  Total  tonnage  to  be  moved  from  base  ‘c’  to  base  ‘cd’  must 
equal  the  direct  deliveries  plus  the  transshipment  deliveries  plus  the  amount 
left  undelivered. 

2.  CARGB(ct,  cd):  The  amount  of  cargo  being  unloaded  at  base  ‘ct’  for  trans¬ 
shipment  to  base  ‘cd’  must  equal  the  amount  being  picked  up  at  ‘ct’  for  trans¬ 
shipment  delivery  to  base  ‘cd’. 
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3.  PLl (routes):  Cargo  carried  on  the  first  leg  of  a  route  cannot  exceed  the  sum 
of  the  payloads  of  all  aircraft  flying  that  route. 

4.  PLP(routes,  stp):  Cargo  loaded  to  fly  leg  ‘stp’  {stp  >  1)  cannot  exceed  the 
unused  capacity  on  the  previous  leg  plus  the  amount  of  cargo  just  unloaded. 

5.  ABAL(c,  aircraft):  At  each  base  ‘c’  the  number  of  landings  by  an  aircraft  type 
must  equal  the  number  of  takeoffs  by  that  type  (conservation  of  flow). 

6.  FREQ(c,  cd):  The  required  minimum  number  of  direct  flights  (i.e.  the  fre¬ 
quency  requirement)  from  base  ‘c’  to  base  ‘cd’  must  be  attained. 

7.  UPHRS (aircraft):  Maximum  number  of  flying  hours  for  each  aircraft  type  must 
not  be  exceeded. 

8.  LOHRS(aircraft):  Minimum  number  of  flying  hours  for  each  aircraft  type  must 
be  achieved. 

9.  MAXLND(c):  The  maximum  allowable  visits  to  a  base  ‘c’  must  not  be  ex¬ 
ceeded. 

The  objective  function,  COST,  is  the  sum  of  (a)  aircraft  operating  hours,  (b)  a 
surcharge  of  79%  on  flying  hours  for  commercial  missions  that  are  bought  one-way, 
(c)  a  handling  charge  for  each  time  a  ton  of  cargo  is  loaded  or  unloaded,  and  a 
“penalty  cost”  of  $10,000  per  ton  for  failure  to  deliver  cargo. 

STORM  Model  in  GAMS 

*  STORM  MODEL 
$OFFSYMXREF  OFFSYMLIST 

OPTIONS  S0LPRINT=0FF,  SYS0UT=0FF,RESLIM=10000  ; 

SETS 

aircraft  PLANE  TYPES  / 

$INCLUDE  "planes. stm" 

/ 
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stp  NUMBERED  SEQUENCE  OF  STOPS  IN  ROUTE  /01  ♦  18/  ; 

ALIAS  (stp,st2); 

SET  routes  AIRCRAFT  ROUTES  / 

$INCLUDE  "route. tmp" 

/  ; 

SET  c  CITIES  / 

$INCLUDE  "bases. stm" 

/  ; 

ALIAS  (c,cd,ct)  ; 

SETS  tp(ct,cd)  TRANSSHIPMENT  BALANCE  PAIRS  / 

$INCLUDE  "trans.tmp" 

/ 

dvar (routes, stp, st2)  COMBINATIONS  FOR  DIRECT  DELIVERY  / 
$INCLUDE  "d.var" 

/ 

svar (routes, stp, st2)  SECOND  HALF  OF  TRANSSHIPMENT  DELIVERY  / 
$INCLUDE  "s.var" 

/ 

tvar (routes, stp, st2,cd)  FIRST  HALF  OF  TRANSSHIPMENT  DELIVERY  / 
$INCLUDE  "t.var" 

/ 

onrte (routes, stp, st2,c,ct)  EACH  ROUTES  SEQUENCE  OF  STOPS  / 
$INCLUDE  "onroute.tmp" 

/ 

start (rout es,c)  STARTING  CITY  OF  EACH  ROUTE  / 

$INCLUDE  "start. tmp" 

/ 

ends (routes, c)  ENDING  CITY  OF  EACH  ROUTE  / 
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$IMCLUDE  "ends.tmp" 

/  ; 

PARAMETERS  LENGTH (routes)  NUMBER  OF  CITIES  VISITED  BY  ROUTE  / 
$INCLUDE  " lengths. tmp" 

/ 

V (routes, c)  NUMBER  OF  VISITS  TO  EACH  CITY  BY  EACH  ROUTE  / 
$INCLUDE  "visits. tmp" 

/ 

TC(c,cd)  TOTAL  CARGO  TO  BE  MOVED  BETWEEN  CITIES  / 
$INCLUDE  "load. tmp" 

/ 

CF(c,cd)  FREQUENCY  REQS  BETWEEN  CITIES  / 

$INCLUDE  "freq.tmp" 

/  ; 


TABLE  H(routes , aircraft)  FLYING  HOUR  ARRAY 
$IWCLUDE  "fy94aug.tms" 


PARAMETER  ML(c)  MAXIMUM  NUMBER  OF  MISSIONS  TO  EACH  CITY  ; 
ML(c)  =1000; 

TABLE  AC (aircraft,*)  AIRCRAFT  INFORMATION 


$INCLUDE  "misc.stm" 
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SCALAR  B  CARGO  HANDLING  COST  PER  TON  lYTQ.Ql  ; 

SCALARS  CARGUNM,  CARGTOT  ; 

VARIABLE  Z  ; 

POSITIVE  VARIABLES 

X (routes, aircraft)  SORTIES  ON  A  ROUTE  BY  AN  AIRCRAFT  TYPE 
Y(c,cd)  UNDELIVERED  CARGO 
D (routes, stp,st2)  DIRECT  DELIVERY  CARGO 

T (routes, stp,st2,cd)  FIRST  PART  OF  TRANSSHIPMENT  DELIVERY 
S (routes, stp,st2)  SECOND  PART  OF  TRANSSHIPMENT  DELIVERY 
SLK(routes,stp)  UNUSED  CAPACITY  ON  EACH  LEG  OF  EACH  ROUTE  ; 

***  BOUNDS  ON  'X'  VARIABLES  INCLUDED  IN  FILE  'LXODA.STM'. 

$INCLUDE  "Ix.stm" 

EQUATIONS 

TREQ(c,cd)  tonnage  required  to  move  between  city  pairs 
CARGB(ct,cd)  cargo  transshipped  in  equals  cargo  transhipped  out 
PLl (routes)  capacity  constraint  on  first  leg  of  a  route 
PLP (routes, stp)  capacity  constraint  on  other  legs  of  route 
ABAL(c, aircraft)  landings  equals  takeoffs  for  each  plane  type 
FREQ(c,cd)  meet  frequency  requirements 
UPHRS (aircraft)  upper  limit  on  aircraft  flying  hours 
L0HRS( aircraft)  lower  limit  on  aircraft  flying  hours 
MAXLND(c)  maximum  number  of  landings  at  any  base 
COSTS  objective  function  of  flying  and  handling  costs  ; 

TREQ(c,cd)$TC(c,cd) . .  Y(c,cd)  + 

SUM( (routes , stp , st2) $ (dvar (routes , stp , st2) $onrte (routes , stp , st2 , c , cd) ) , 
D(routes,stp,st2))  + 


A-7 


SUM ( (routes , stp , st2 , ct ) $ (tvar (routes , stp , st2 , cd) $onrte (routes , stp , st2 , c , ct ) ) , 
T (routes, stp,st2,cd))  =G=  TC(c,cd)  ; 

CARGB(ct,cd)  $tp(ct,cd).. 

SUM ( (routes , stp , st2 , c)  $ (tvar (routes , stp , st 2 , cd) $onrte (routes , stp , st2 , c , ct ) ) , 
T (routes, stp, st2,cd))  =E= 

SUM( (routes , stp , st2)  $ (svar (routes , stp , st2) $onrte (routes , stp , st2 , ct , cd) ) , 

S (routes, stp, st2))  ; 

PLl(routes) . .  SUM(aircraft  $H(routes .aircraft) , 

AC (aircraft, "PAY")  *  X (routes , aircraft) )  =E=  SLK(routes, "01")  + 

SUM(st2  $dvar(routes , "01" , st2) ,  D(routes,"01" ,st2))  + 

SUM(st2  $svar (routes, "01", st2) ,  S (routes, "01", st2))  + 

SUM((st2,cd)  $tvar(routes,"01",st2,cd) ,  T(routes,"01" ,st2,cd))  ; 

PLP (routes, stp) $(0RD (stp)  It  LENGTH (routes)  and  ORD(stp)  ge  2).. 
SLK(routes,stp-l)  + 

SUM(st2  $dvar (routes, st2, stp) ,  D (routes, st2, stp))  + 

SUM(st2  $svar (routes, st2, stp) ,  S (routes, st2, stp))  + 

SUM((st2,cd)  $tvar (routes, st2, stp, cd) ,  T (routes, st2, stp, cd)) 

=E=  SLK (routes , stp)  + 

SUM(st2  $dvar (routes , stp, st2) ,  D(routes,stp,st2))  + 

SUM(st2  $svar (routes, stp, st2) ,  S (routes, stp, st2))  + 

SUM((st2,cd)  $tvar (routes, stp, st2, cd) ,  T (routes, stp, st2,cd))  ; 

ARAL (c, aircraft) . . 

SUM(routes$ (start (routes , c)  $H(routes , aircraft) ) , X(routes , aircraft) ) 
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=E= 

SUM (rout es$ (ends (routes, c)  $H(routes, aircraft)) ,X(routes, aircraft))  ; 

FREQ (c , cd)  $CF (c , cd) . . 

SUM( (routes .aircraft )$H (routes, aircraft) , 

SUM ( ( stp , st2) $onrt e (rout es , stp , st2 , c , cd) , 

X(routes, aircraft)))  =G=  CF(c,cd)  ; 

UPHRS (aircraft) . .  SUM(routes  $H(routes,aircraft) , 

H(routes .aircraft)  *  X(routes,aircraft)) 

=L=  AC(aircraft,"FH")  ; 

LQHRS (aircraft) . .  SUM(routes  $H(routes,aircraft) , 

H(routes, aircraft)  *  X(routes, aircraft)) 

=G=  AC(aircraft,"LH")  ; 

MAXLND(c)$(ML(c)  le  1000)..  SUM ((routes, aircraft) 

$H(routes, aircraft) ,  V(routes,c)  *  X(routes .aircraft) )  =L=  ML(c)  ; 

COSTS..  Z  =E=  SUM( (routes, aircraft)  $H(routes, aircraft) , 

AC (aircraft, "COST")  *  H (routes, aircraft)  *  X (routes .aircraft) ) 

+  SUM( (routes, aircraft)  $(start(routes,"EXXX")  $H(routes,aircraft)) , 

0.79  *  AC (aircraft , "COST")  *  H(routes, aircraft)  *  X (routes, aircraft)) 

+  B  *  SUM((routes,stp,st2)$dvar(routes,stp,st2) ,D(routes,stp,st2)) 

+  2  *  B  *  SUM( (routes, stp, st2)$svar (routes, stp, st2) ,S (routes, stp, st2)) 
+  20000  *  SUM((c,cd)$TC(c,cd),Y(c,cd))  ; 

MODEL  STORMl  /ALL/  ; 
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OPTION  ITERLIM  =  40000  ; 

OPTION  LP  =  ZOOM; 

SOLVE  STORMl  USING  LP  MINIMIZING  Z  ; 
DISPLAY  Z.L  ; 
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Appendix  B.  TREQ  and  FREQ  Input  Files  For  STORM 

TREQ  Data  File _ 


KDOV.LETO 

22.12 

KDOV.LIPA 

173.77 

KSUU.FJDG 

40.33 

LETO.KDQV 

58.26 

EDAR.KWRI 

9.60 

EDAR.KNGU 

23.60 

EDAR.KDOV 

947.49 

EDAF.KDOV 

923.74 

KDOV.EDAR 

744 . 00 

KDOV.EDAF 

787 . 57 

KDOV.LTAG 

140.19 

RJSM.KSUU 

57.58 

KSUU.RJSM 

108.69 

PHIK.KSUU 

451.39 

KSUU.PHIK 

545.56 

KDOV.OERY 

199.12 

KDOV.LGIR 

9.07 

KSUU.RKJK 

32.71 

KSUU.RKPK 

20.47 

EGUN.KCHS 

345.78 

OEDR.KDOV 

106.79 

KDOV . OEDR 

323.47 

PKWA.KSUU 

27.65 

KSUU.RODN 

392.99 

KSUU.RKSO 

634.20 
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KSUU.PGUA 

196.14 

KSUU.RJTY 

340 . 32 

RODN.KSUU 

254.33 

RJTY.KSUU 

338.31 

KSUU.PAED 

239.52 

KSUU.RJOI 

25.02 

KSUU.PJON 

11.91 

KTCM.PAED 

62.50 

KTCM.PADK 

81.12 

PADK.KTCM 

9.96 

RJTY . PHIK 

13.83 

RJTY.PGUA 

25.51 

RJTY.WSAP 

28.43 

KSUU.WSAP 

18.68 

WSAP.KSUU 

6.92 

RJTY.FJDG 

162.67 

RJTY.RKJK 

0.49 

RJTY.RKSO 

50.90 

RJTY.RJFF 

13.58 

RJTY.RJOI 

27.51 

RJTY.RJAW 

8.71 

RJTY.RODN 

69.03 

RJTY.RJSM 

159.37 

RJTY.RJCB 

1.66 

RJTY.RJAM 

2.68 

PAED.KSUU 

330.32 

PHIK.PGUA 

17.56 

PHIK. ROOM 

17.69 

PHIK.PJON 

163.57 

PHIK.PKWA 

122.46 

PHIK.PWAK 

70.38 

PHIK.PMDY 

11.82 

PGUA.KSUU 

83.91 

NZCH.KSUU 

20.14 

FJDG.KSUU 

95.26 

FJDG.RJTY 

27.14 

FJDG.OMFJ 

15.76 

OMFJ.FJDG 

0.01 

FJDG.WSAP 

17.10 

WSAP.FJDG 

12.69 

WSAP.RJTY 

2.31 

RKSO.KSUU 

647.15 

RKSO.RJTY 

28.92 

RKSO.RODN 

25.70 

RJFF.RJTY 

0.07 

RJOI.KSUU 

14.30 

RJOI.RODN 

26.18 

RJAW.RJTY 

11.86 

RODN.RJTY 

87.59 

RODN.PHIK 

25.37 

RODN.RKPK 

2.91 

RODN.RKJK 

3.69 

RODN.RKSO 

41.46 

RODN.RJOI 

30.77 

RODN.RJSM 

3.68 

RODN.PGUA 

20.17 

PGUA.RODN 

9.89 

RJSM.RJTY 

52.20 

RJSM.RKJK 

0.03 

RJSM.RODN 

2.10 

RJCB.RJTY 

1.77 

ABAS.KSUU 

18.58 

APLM.ASRI 

0.01 

APLM.KSUU 

31.68 

ASRI.APLM 

2.57 

ASRI.KSUU 

50.51 

APWR.KSUU 

17.24 

RJAM.RJTY 

1.00 

KSUU.NZCH 

22.55 

KSUU.ABAS 

106.86 

KSUU.ASRI 

34.67 

KSUU.APWR 

11.04 

PJON.PHIK 

30.43 

PKWA.PHIK 

73.30 

PWAK . PHIK 

0.39 

PMDY.PHIK 

4.49 

KSUU.VTBD 

9.09 

RODN.VVMB 

0.01 

VVNB.RODN 

0.01 

RODN.VTBD 

17.53 

KCHS . SBGL 

1.59 

SBGL.KCHS 

6.94 

KCHS . SEQU 

14.81 

RJFF.KSUU 

6.41 

KSUU.RJFF 

16.32 

KDOV.RODN 

8.87 

KCHS.SVMI 

9.31 

SKBO.KCHS 

12.01 

KSUU.WIIH 

5.06 

WIIH.KSUU 

4.33 

KCHS.MHTG 

17.25 

MHTG.KCHS 

6.68 

KCHS.EGUN 

292.91 

KCHS.MHSC 

216.94 

KCHS.MPHQ 

493.74 

KCHS.MSSS 

52.38 

KDOV.OJAF 

8.64 

KDQV.HECA 

55.52 

KNGU.LERT 

132.13 

KNGU.LIRN 

64.63 

KNGU.LICZ 

105.19 

KNGU . OBBI 

63.69 

KNGU.BIKF 

142.76 

KNGU.MUGM 

271.30 

KNGU.TJNR 

211.93 

KWGU.TXKF 

128.44 

KNGU.OMFJ 

40.88 

OMFJ.KNGU 

0.01 

KWRI.LPLA 

285.47 

KWRI.BGTL 

146.24 

EDAR.EGUM 

58.06 

EDAR.LETO 

7.40 

EDAR.LTAG 

17.96 

EDAR.HECA 

21.67 

EDAR.LIPA 

51.65 

EDAF.EGUN 

81.32 

EDAF.LETO 

5.97 

EDAF.OEDR 

162.35 

EDAF.LTAG 

211.83 

EGUN.EDAR 

27.97 

EGUN.LETO 

3.74 

EGUN.LIPA 

4.84 

EGUN.LTAG 

18.42 

LERT.KNGU 

138.37 

LERT.LICZ 

12.32 

LETO.EDAR 

2.31 

LETO.EDAF 

10.02 

LETO.EGUN 

0.88 

LETO.LIPA 

16.57 

LPLA.KWRI 

99.07 

LIPA.EDAR 

0.05 

LIPA.EGUN 

0.55 

LIPA.LETO 

0.02 

LIRN.KNGU 

0.01 

LIRN.LIEO 

0.01 

KNGU.LIPA 

21.26 

LIRN.LGSA 

0.01 

LGSA.LIRN 

0.01 

LIEO.LIRN 

0.01 

LIRN.LICZ 

0.01 

LICZ.KNGU 

143.87 

LICZ.LIRN 

18.47 

LTAG.KDOV 

176.99 

LTAG.EDAR 

101.43 

LTAG.EDAF 

60.28 

LTAG.EGUN 

71.51 

OJAF.KDOV 

5.19 

OBBI.KNGU 

110.03 

OEDR.EDAF 

15.09 

BGTL.KWRI 

156.18 

MHSC.KCHS 

83.21 

MHSC.MPHO 

29.94 

MPHQ.KCHS 

356.59 

MPHO.MHSC 

109.61 

MUGM.KNGU 

34.81 

TJNR.KNGU 

42.64 

TJNR.MTPP 

0.01 

TJNR.TAPA 

1.07 

TAPA.TJNR 

0.01 

TJNR.TISX 

0.01 

TISX.TJNR 

0.01 

TXKF.KNGU 

28.11 

KWRI.CYYT 

0.01 

CYYT.KWRI 

0.01 

KCOF.FHAW 

77.06 

KCOF.TAPA 

52.91 

TAPA.KCOF 

6.53 

KNGU.HKNA 

3.70 

PAED.PADK 

10.15 

PHIK.ASRI 

13.12 

KCHS . SLLP 

37.02 

MPHO . SGAS 

1.40 

SGAS.MPHQ 

0.16 

MPHO . SKBO 

9.42 

SKBO.MPHO 

0.35 

MPHO . SAEZ 

4.59 

SAEZ.MPHQ 

0.04 

MPHO.MGGT 

8.39 

MGGT.MPHQ 

0.02 

MPHO . SPIM 

12.90 

SPIM.MPHO 

0.60 

MPHO . SLLP 

14.02 

SLLP. MPHO 

0.89 

MPHO.MNMG 

14.06 

MPHO , SVMI 

0.71 

MPHO . SUMU 

0.03 

MPHO.MROC 

3.20 

MPHO . SBGL 

0.74 

MPHO. MS SS 

23.18 

MSSS.MPHO 

0.01 

MPHO.SCEL 

0.65 

SCEL.MPHO 

0.01 

MPHO.MHTG 

19.69 

MHTG.MPHO 

0.13 

MPHO . SEQU 

12.63 

SEQU.MPHO 

0.50 

RJTY.PAED 

16.60 

MPHO.SBBR 

6.22 

PWAK.KSUU 

6.53 

RKJK.KSUU 

25.67 

HECA.KDOV 

16.39 

HECA.EDAR 

5.62 

KCHS.SGAS 

1.17 

KCHS.SAEZ 

0.10 

SAEZ.KCHS 

0.01 

KCHS.GOOY 

10.64 

KCHS.FZAA 

2.55 

SLLP.KCHS 

0.39 

KCHS . SUMU 

0.32 

SUMU.KCHS 

0.01 

KCHS. FIT J 

4.24 

KCHS.SCEL 

9.18 

SCEL.KCHS 

11.92 

LIPA.KDOV 

3.41 

KDOV.LLBG 

5.65 

PAED.RJTY 

18.74 

MROC.MPHO 

3.00 

SBGL.MPHO 

7.01 

RKSO.RKJK 

10.88 

EDAR.OJAF 

19.14 

EDAR.OEDR 

14.08 

KSUU.NSTU 

1.10 

PJON.KSUU 

18.63 

KSUU.PKWA 

70.41 

PAED.KTCM 

14.41 

MNMG.MPHQ 

0.01 

SVMI.MPHO 

0.23 

SUMU.MPHO 

0.09 

KCHS.GLRB 

0.34 

GLRB.KCHS 

0.01 

KCHS.MUGM 

14.82 

KCHS.DRRN 

5.53 

KCHS.TJNR 

8.68 

KCHS.SPIM 

14.93 

SPIM.KCHS 

3.67 

KDOV . EGUN 

6.08 

KTCM.RQDN 

10.71 

KTCM.PHIK 

6.39 

SBBR.KCHS 

7.13 

PHIK.KTCM 

51.49 

EGUN. KDOV 

35.71 

RJTY.KTCM 

31.05 

RKSO.KTCM 

14.96 

LICZ.KDOV 

9.81 

EDAF.LIPA 

69.31 

EDAF . OERY 

100.03 

KCHS . SKBO 

64.34 

KCHS.SBBR 

26.89 

KCHS. EDAF 

3.91 

KCHS.MNMG 

0.22 

EDAF.LIRN 

1.84 

EDAF.LGIR 

10.46 

KNGU.TISX 

29.87 

KSUU.PWAK 

1.87 

KSUU.OBBI 

0.04 

KSUU.OEDR 

0.02 

KTCM.RJTY 

0.12 

KWRI.EDAR 

1.40 

FHAW.KCOF 

18.09 

EGUN.KNGU 

0.19 

TJNR.KCHS 

36.62 

LTAG.LIPA 

7.65 

LTAG.LIRN 

6.23 

LTAG.LLBG 

1.31 

LTAG.LETG 

2.13 

PADK.PAED 

2.33 

PADK.PHIK 

1.92 

LIPA.EDAF 

0.03 

VTBD.RODN 

2.76 

WIIH.RODN 

0.01 

RODN.WIIH 

1.96 

RODN.FJDG 

2.67 

RODN.WSAP 

10.55 

PAED.PHIK 

4.60 

PHIK.PAED 

2.95 

PHIK.RJTY 

7.57 

PHIK.RKSO 

3.02 

PHIK.NSTU 

1.84 

RJOI.RJTY 

6.99 

RKJK.RODN 

8.22 

RKJK.RJTY 

1.80 

RKJK.RKSO 

3.81 

APLM.PHIK 

0.35 

RJSM.PHIK 

0.19 

FJDG.RODN 

0.47 

TJNR.TBPB 

8.54 

TJMR.MUGM 

0.24 

RJTY.VTBD 

2.57 

RJTY.EDAF 

0.28 

RJTY.RKPK 

0.16 

RJTY.RKTN 

0.18 

MKJP.MUGM 

6.21 

LIPA.LIRN 

4.09 

RKSO.PHIK 

4.19 

RKSO.RKPK 

12.80 

NSTU.PHIK 

0.40 

ASRI.ABAS 

2.70 

ASRI.APWR 

6.88 

EDAR.LPLA 

20.44 

EGUN . EDAF 

18.93 

EGUN.LPLA 

0.09 

BIKF.KNGU 

2.86 

EDAR.LIRN 

6.83 

EDAR.LICZ 

16.59 

EDAR.LGIR 

1.34 

LETO.LTAG 

2.71 

LETO.OBBI 

0.23 

LETO.LIRN 

3.09 

LETO.LERT 

5.60 

LETO.LICZ 

1.95 

LETO.LGIR 

1.99 

LERT.LPLA 

0.42 

LERT.LIRN 

8.19 

LICZ.OBBI 

11.51 

LICZ.OMFJ 

1.22 

LICZ . OEDR 

1.98 

LICZ.LETQ 

0.07 

LICZ.LERT 

11.36 

LICZ.LTAG 

4.82 

PGUA.PHIK 

4.93 

PGUA.RJSM 

0.01 

PGUA.RJTY 

2.83 

PGUA.RKSQ 

9.66 

LETO.LLBG 

9.90 

LERT.LTAG 

5.51 

EGUN.LICZ 

5.85 

EDAF.OJAF 

5.00 

KCHS.MGGT 

17.36 

LICZ.LIPA 

11.99 

LGIR.KDOV 

19.61 

KCHS.TAPA 

18.67 

KSUU.APLM 

25.32 

PADK.KSUU 

29.37 

SGAS.KCHS 

7.16 

PHIK.FJDG 

9.79 

SVMI.KCHS 

5.55 
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APLM.PGUA 


5.74 


RJTY.PWAK 

4.05 

MSSS.KCHS 

27.50 

LLBG.KDOV 

10.28 

LTAG.LGIR 

4.02 

PAHT.PAED 

4.53 

OJAF.EDAR 

6.34 

OBBI.LICZ 

9.34 

RODN.PAED 

6.68 

PAED.RODN 

7.99 

PAED.RKSO 

25.54 

EDAF.EDAR 

4.64 

EDAF.LICZ 

8.16 

EDAF.LLBG 

11.01 

MPHO.MZBZ 

9.79 

RJTY.OMFJ 

45.88 

ASRI.PHIK 

5.01 

EDAR.EDAF 

4.86 

LERT . OBBI 

4.52 

LICZ.EDAF 

4.03 

LGIR.LIPA 


7.53 


FREQ  Data  File 


KSUU.PGUA 

4.29 

KTCM.PADK 

4.29 

PADK.KTCM 

4.29 

RJTY.WSAP 

16.5 

RJTY.FJDG 

4.29 

RJTY.RKJK 

4.29 

RJTY.RKSO 

13.5 

RJTY.RJFF 

4.29 

RJTY.RJAW 

4.29 

RJTY.RJSM 

13.5 

RJTY.RJCB 

4.29 

RJTY.RJAM 

4.29 

PHIK.PJON 

8.57 

PHIK.PKWA 

8.57 

PHIK.PWAK 

6.43 

PHIK . PMDY 

4.29 

PGUA.KSUU 

4.29 

FJDG.RJTY 

4.29 

F JDG . OMFJ 

12.86 

OMFJ.FJDG 

12.86 

FJDG.WSAP 

13.5 

WSAP.FJDG 

13.5 

WSAP.RJTY 

16.5 

RKSO.RJTY 

13.5 

RKSQ.RODN 

13.5 

RJFF.RJTY 

4.29 

RJAW.RJTY 

4.29 

RODN.RKPK 

4.29 
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RODN.RKJK 

13.5 

RODN.RKSO 

13.5 

RODN.RJSM 

13.5 

RODN.PGUA 

8.57 

PGUA.RODN 

8.57 

RJSM.RJTY 

13.5 

RJSM.RKJK 

13.5 

RJSM.RODN 

13.5 

RJCB.RJTY 

4.29 

APLM.ASRI 

4.29 

ASRI.APLM 

4.29 

RJAM.RJTY 

4.29 

PJON.PHIK 

8.57 

PKWA.PHIK 

8.57 

PWAK.PHIK 

6.43 

PMDY.PHIK 

4.29 

RODN.VVNB 

1.07 

VVNB.RODW 

1.07 

KCHS . SBGL 

2.14 

SBGL.KCHS 

2.14 

KCHS.MHSC 

4.29 

KCHS.MPHO 

4.29 

KDOV.OJAF 

2.14 

KNGU.LIRN 

12.86 

KNGU.LICZ 

15.5 

KNGU.OBBI 

12.86 

KNGU.MUGM 

12.86 

KMGU.TJNR 

4.29 

KNGU.TXKF 

3.21 

KNGU.OMFJ 

15.5 

OMFJ.KNGU 

15.5 

KWRI.LPLA 

13.5 

KWRI.BGTL 

4.29 

EDAR.EGUN 

12.86 

EDAR.LETO 

4.29 

EDAR.HECA 

4.29 

EDAR.LIPA 

8.57 

EDAF.LETQ 

4.29 

EDAF . OEDR 

24.5 

EDAF.LTAG 

12.86 

EGUN . EDAR 

12.86 

EGUN.LETO 

4.29 

EGUN. LI PA 

4.29 

LETO . EDAR 

4.29 

LETO . EDAF 

4.29 

LETO . EGUN 

4.29 

LETO. LI PA 

4.29 

LPLA.KWRI 

13.5 

LI PA. EDAR 

8.57 

LIPA.EGUN 

4.29 

LIPA.LETO 

4.29 

LIRN.KNGU 

12.86 

LIRN.LIEO 

12.86 

LIRN.LGSA 

4.29 

LGSA.LIRN 

4.29 

LIEO.LIRN 

12.86 

LIRN.LICZ 

12.86 

LICZ.KNGU 

15.5 

LICZ.LIRN 

12.86 

LTAG.EDAF 

12.86 

OJAF.KDOV 

2.14 

OBBI.KNGU 

12.86 

□EDR.EDAF 

24.5 

BGTL.KWRI 

4.29 

MHSC.KCHS 

4.29 

MHSC.MPHO 

8.57 

MPHO.KCHS 

4.29 

MPHO.MHSC 

8.57 

MUGM.KNGU 

12.86 

TJNR.KNGU 

4.29 

TJNR.MTPP 

1.07 

TJNR.TAPA 

2.14 

TAPA.TJNR 

2.14 

TJNR.TISX 

4.29 

TISX.TJNR 

4.29 

TXKF.KNGU 

3.21 

KWRI . CYYT 

2.14 

CYYT.KWRI 

2.14 

KCHS . SLLP 

2.14 

MPHO . SGAS 

2.14 

SGAS.MPHO 

2.14 

MPHO . SKBO 

1.07 

SKBO.MPHO 

1.07 

MPHO.SAEZ 

2.14 

SAEZ.MPHO 

2.14 

MPHO.MGGT 

2.14 

MGGT.MPHO 

2.14 

MPHO.SPIM 

2.14 

SPIM.MPHO 

2.14 

MPHO . SLLP 

2.14 

SLLP.MPHQ 

2.14 

MPHO.MNMG 

2.14 

MPHO . SVMI 

1.07 

MPHO . SUMU 

2.14 

MPHO.MROC 

2.14 

MPHO . SBGL 

2.14 

MPHO.MSSS 

8.57 

MSSS.MPHO 

8.57 

MPHO . SCEL 

2.14 

SCEL.MPHO 

2.14 

MPHO.MHTG 

8.57 

MHTG.MPHO 

8.57 

MPHO.SEQU 

1.07 

SEQU.MPHO 

1.07 

MPHO . SBBR 

2.14 

HECA . EDAR 

4.29 

KCHS.SGAS 

2.14 

KCHS . SAEZ 

2.14 

SAEZ.KCHS 

2.14 

KCHS.FZAA 

2.14 

SLLP. KCHS 

2.14 

KCHS . SUMU 

2.14 
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SUMU.KCHS 


2.14 


KCHS.FTTJ 

1.07 

KCHS.SCEL 

2.14 

SCEL.KCHS 

2.14 

MROC.MPHO 

2.14 

SBGL.MPHO 

2.14 

MNMG.MPHO 

2.14 

SVMI.MPHQ 

1.07 

SUMU.MPHO 

2.14 

KCHS.GLRB 

2.14 

GLRB.KCHS 


2.14 


Appendix  C.  Application  of  Simple  Point  Kriging  to  the  Residuals  of 

Metamodel  1 

In  kriging,  the  first  step  is  to  determine  the  distance  matrix,  h.  Here,  the 
distance  matrix  was  obtained  calculating  the  Euclidean  distances  between  the  design 
points.  The  Euclidean  distance  is  given  by 

hi  =  \J{xi  -  +  (yi  -  yof  (C.l) 

where  x  and  y  are  the  coordinate  values  of  the  design  points.  Recalling  the  first 
order  2^  full  factorial  design  of  Metamodel  1  in  Chapter  4,  the  following  distance 
matrix  is  obtained; 

022  2^/2 

2  0  2v^  2 

h  = 

2  2^2  0  2 
2^/2  2  2  0 

Also,  we  can  derive  the  distances  of  the  points  that  are  to  be  estimated  to  the 
other  points.  For  example,  the  distance  vector  of  the  first  validation  point  in  the 
validation  design  (-0.5,  -0.5)  to  the  test  design  points  is 

0.707 

1.581 

1.581 

2.121 

The  second  step  is  to  construct  the  Equation  (2.9)  using  a  semi-veriogram 
function  model.  Here,  the  spherical  model  given  in  Equation  (2.10)  was  used.  The 
constant  parameters  are  choosen  as  Co  =  0,  C  =  1.444  x  10^*^  which  is  the  vari¬ 
ance  between  regression  metamodel  residuals,  and  a  =  2\/2  which  is  the  maximum 
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distance  in  the  distance  matrix.  By  substituting  the  distance  values  into  spherical 
model,  for  the  fist  validation  point,  Equation  (2.9)  becomes 


0 

1 

1 

1 

0 

1.011x101° 

1 

1.011x101° 

0 

1 

1.011x101° 

1.144x101° 

1 

1.144x101° 

1.011x101° 

1 

1 

1.011x101° 

1.144x101° 

1.144x101° 

1.011x101° 

0 

1.011x101° 

1.011x101° 

0 

A 

1 

u;i 

4.201x10° 

102 

= 

8.594x10° 

LO3 

8.594x10° 

LO4 

1.046x101° 

The  Lagrangian  multiplier  and  the  weights  are  calculated  as 

4.725x10^ 

0.596 

0.178 

0.178 

0.048 

Finally,  substituting  the  weights  into  Equation  (3.12),  we  find  the  kriged  resid¬ 
ual  for  the  first  validation  design  point  as  Ci  =  —26680.  The  rest  of  the  kriged 
residuals  were  found  by  following  the  same  steps. 
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Appendix  D.  The  Channel  Cargo  System  of  Air  Mobility  Command 


Figure  D.l  World  Map  Showing  Cargo  Channels 
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